The aim of this website is to report the ongoing geospatial data analysis for the CCRI’s Agricultural Transition: Heart of South West LEP project. All of the tables and plots shown on this website can also be accessed as individual files HERE. The full R programming code used to perform the work is included in this report, and can be accessed by clicking on the “Code” buttons. If you are unfamiliar with R code, comments (in plain English) have been added to describe the operation of each code block, so it should be possible to follow the methodology.


1. Pillar 1: County-level

This section shows the code and outputs for the analysis of Pillar 1 (P1) payments at county level (for the counties of Cornwall, Devon, Dorset, and Somerset)


1.1 Data import & cleaning

Import and cleaning of the CAP Payments (data source: HERE)

Step 1 - Load code packages and import data

#> 1.1.1 Load libraries 
library(tidyverse) # data munging and analysis
library(sf) # simple features for GIS
library(qgisprocess) # access to QGIS algorithms
library(here) # relative path management for reproducibility
library(lubridate) # for date string manipulation/conversion
library(janitor) # data cleaning functions
library(knitr) # report rendering with rmarkdown
library(ggplot2) # plots and visualisations
library(readxl) # import xls
library(scales) # for plot scale customisation
library(mapproj) # map zoom coords
library(leaflet) # interactive web maps
library(reshape2) # to reshape data (used for "melt" function)
  
#> 1.1.2 Import data (trimming white space)
#> "RPA" sheet
rpa.1 <- readxl::read_xlsx(here("In", "BPS_Glos", "2020_All_CAP_Search_Results_Data_P14.xlsx"), trim_ws = TRUE, sheet = "RPA")
#> "RPA2" sheet
rpa.2 <- readxl::read_xlsx(here("In", "BPS_Glos", "2020_All_CAP_Search_Results_Data_P14.xlsx"), trim_ws = TRUE, sheet = "RPA2")

Step 2 - Data cleaning and pre-processing

#> 1.1.3 Data cleaning and pre-processing

#> Bind the two RPA tables together (they are split alphabetically by beneficiary name)
db <- rbind(rpa.1, rpa.2)

#> Remove duplicates
db <- db %>% 
  distinct()

#> Remove any empty records where all columns are "NA"
# Following filters rows with at least one column not "NA"
db <- janitor::remove_empty(db, which = "rows")

#> Select only relevant rcolumns
db <- db %>% 
  select("BeneficiaryCode", "PostcodePrefix_F202B", "TownCity_F202C", "Basic payment scheme", "Greening: practices beneficial for climate and environment") |> 
  arrange(desc(`Basic payment scheme`))

Table 1.1.1 Full database on import = 88443 records (all England)

## Rows: 88,443
## Columns: 5
## $ BeneficiaryCode                                              <dbl> NA, NA, N~
## $ PostcodePrefix_F202B                                         <chr> "LN4", "S~
## $ TownCity_F202C                                               <chr> "LINCOLN"~
## $ `Basic payment scheme`                                       <dbl> 1929366.7~
## $ `Greening: practices beneficial for climate and environment` <dbl> 872521.2,~

Table 1.1.2 First 6 rows of database

## # A tibble: 6 x 5
##   BeneficiaryCode PostcodePrefix_F202B TownCity_F202C `Basic payment scheme`
##             <dbl> <chr>                <chr>                           <dbl>
## 1              NA LN4                  LINCOLN                      1929367.
## 2              NA SN2                  SWINDON                      1838669.
## 3              NA CB2                  Cambridge                     975272.
## 4              NA NE71                 WOOLER                        969240.
## 5              NA LN5                  LINCOLN                       875713.
## 6              NA PE38                 DOWNHAM MARKET                847869.
## # ... with 1 more variable:
## #   `Greening: practices beneficial for climate and environment` <dbl>


1.2 Land cover processing

The CAP Payments dataset reports payments at postcode district level, but postcode districts do not conform to administrative boundaries (e.g. counties, district unitary authorities). Our approach is to calculate the area of agricultural land (using Corine land cover data) within each postcode district, and to use this to calculate the P1 payments. For districts which straddle county boundaries, the area of agricultural land both inside and outside of the target county are calculated, and the P1 payments are calculated based on the proportion of agricultural land inside the boundary. Agricultural land cover classes were extracted from the Corine 2018 land cover dataset and used for the agricultural area calculations.

The following geospatial datasets are used for this analysis:
- County boundaries - (OS BoundaryLine)
- Postcode district areas - (OpenDoor Postcode Districts)
- Land Cover data - (Corine 2018)

#> Import county boundaries
#> Import ceremonial counties polygon layer from OS BoundaryLine dataset
counties.sw <- st_read(here("In", "Shape", "South_West.UPDATED.2.shp"), quiet = TRUE)

# glimpse(counties.sw)

#> Remove the deleted Bath and NE Somerset polygon (result of pre-geoprocessing in QGIS)
counties.sw <- counties.sw |> 
  filter(is.na(NAME_2))

#> Import postcode district boundary polygon data (source: https://www.opendoorlogistics.com/downloads/)
pcodes.gb <- st_read(here("In", "Shape", "Pcode_Districts_OpenDoor_2017.shp"), quiet = TRUE)

#> Import CORINE 2018 land cover data (polygons)
#> Data source: https://catalogue.ceh.ac.uk/documents/084e0bc6-e67f-4dad-9de6-0c698f60e34d
corine.gb <- st_read(here("In", "Shape", "corine_2018_GB.shp"), quiet = TRUE)
#> Select only agri land classes:
# 211 - Non-irrigated arable land
# 212 - Permanently irrigated land
# 213 - Rice fields
# 221 - Vineyards
# 222 - Fruit trees and berry plantations
# 223 - Olive groves
# 231 - Pastures
# 241 - Annual crops associated with permanent crops
# 242 - Complex cultivation patterns
# 243 - Land principally occupied by agriculture with significant areas of natural vegetation
# 244 - Agro-forestry areas
# 321 - Natural grasslands
# 322 - Moors and heathland
# 412 - Peatbogs
agri.land.classes <- c("211", "212", "213", "221", "222", "223", "231", "241", "242", "243", "244", "321", "322", "412")
#> Subset data based on land cover codes
corine.agri <- corine.gb |> 
  select(ID, CODE_18) |> 
  filter(CODE_18 %in% agri.land.classes)




  
#> Export agri data for detailed visual checking in QGIS
#> Remove previous shapefile export
unlink(here("Out", "Tests", "Corine_Agri.shp"))
unlink(here("Out", "Tests", "Corine_Agri.dbf"))
unlink(here("Out", "Tests", "Corine_Agri.prj"))
unlink(here("Out", "Tests", "Corine_Agri.shx"))
st_write(corine.agri, here("Out", "Tests", "Corine_Agri.shp"), quiet = TRUE)

#> Show an example plot of agri land extracted from Corine
#> Get Cornwall boundaries
cornwall.sf <- counties.sw |>
  filter(NAME == "Cornwall")
#> Clip corine.gb to extent of Cornwall
corine.clip <- st_intersection(corine.agri,cornwall.sf)
# 
# #> Previous static map (superseded by Leaflet map)
# #> Plot
# plot0 <- ggplot(data = cornwall.sf) +
#   ggtitle("Agricultural land: Cornwall (extracted from Corine 2018)") +
#   # theme_bw() +
#   geom_sf() +
#   geom_sf(data = corine.clip, fill = "green")
# #> Render on web page
# print(plot0)

Figure 1.2.1 Extent of agricultural land (derived from Corine 2018)

#> This code chunk produces an interactive leafelt map of the SW region with Corine agri land shown

#> Subset counties.gb, selecting only Cornwall (inc Isles of Scilly), Devon, Somerset, and Dorset
#> Will also use this string in Section 1.3 loop
#> Create filter string (get names first!)
counties <- c("Cornwall", "Devon", "Dorset", "Somerset")
#> Create sf object of 4 counties
south.west <- counties.sw |>
  filter(NAME %in% counties)
#> Import clipped (to SW region) dissolved layer of Corine agri
#> Pre-processed in QGIS and imported as process in R takes too long
corine.sw <- st_read(here("In", "Shape", "Corine.sw.clip.diss.simp.10.diss2.shp"), quiet = TRUE)
#> Convert layers to WGS84 for use with Leaflet
south.west <- st_transform(south.west, 4326)
corine.sw <- st_transform(corine.sw, 4326)
#> Create Leaflet map
leaflet() |>
  addPolygons(data = corine.sw, fillOpacity = 0.5, color = "green", stroke = FALSE) |>
  addPolygons(data = south.west, fillOpacity = 0, label = TRUE) |>
  addTiles()

1.3 P1 Reduction calculations

With all the required input data in place, P1 payments at county level can now be computed. A brief overview of the workflow:

  • Step 1: Select the four target counties from the GB county data layer
  • Step 2: Initiate a programming loop based on the counties (i.e. same code block is run for each county in turn)
  • Step 3: Create a spatial layer of postcode districts that spatially intersect with the “active” county
  • Step 4: Create a spatial layer of postcode districts clipped exactly to the extent of the “active” county
  • Step 5: Create a spatial layer of postcode districts clipped exactly to the extent of the “active” county
  • Step 6: Using the Corine data on agricultural land area (prepared in Section 1.2), calculate the agricultural land area totals for intersected postcode district polygons, and clipped postcode district polygons respectively. This gives us two values: 1) the total area of agricultural land in each postcode district, and; 2) the area of agricultural land in each postcode district that is within the active county
  • Step 7: Using the agricultural land area proportions and the P1 reduction figures, calculate the “baseline” value for P1 in 2020 and following years with payments applied up to and including 2027
  • Step 8: Display output tables and bar plots
  • Step 9: Export output tables and bar plots to CCRI shared drive
#> Create filter string (get names first!)
counties <- c("Cornwall", "Devon", "Dorset", "Somerset")

#> Sequential integer counter (for figure numbers)
i <-0 

#> Initiate loop
for(active_county in counties){
  
#> Loop counter (for table caption number)
i <- i+1

#> Select target county (this will be start of loop)
county <- counties.sw |> 
  filter(NAME == active_county)

#> Select postcode district areas that intersect with active county polygon
pcodes.int <- pcodes.gb[county,]
#> Add area (km2) column
pcodes.int$AreaM2 <- st_area(pcodes.int)
#> Convert from m2 to km2
pcodes.int$AREA_KM2_TOTAL <- pcodes.int$AreaM2 / 1000000
pcodes.int$AREA_KM2_TOTAL <- round(pcodes.int$AREA_KM2_TOTAL, digits = 3)
pcodes.int$AREA_KM2_TOTAL <- as.numeric(pcodes.int$AREA_KM2_TOTAL)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "pcodes.int.shp"))
unlink(here("Out", "Tests", "pcodes.int.dbf"))
unlink(here("Out", "Tests", "pcodes.int.prj"))
unlink(here("Out", "Tests", "pcodes.int.shx"))
st_write(pcodes.int, here("Out", "Tests", "pcodes.int.shp"), quiet = TRUE)


#> Generate a layer of clipped (cookie cutter) postcodes using active county boundary
pcodes.clip <- st_intersection(pcodes.int, county)
#> Add area (km2) column
pcodes.clip$AreaM2 <- st_area(pcodes.clip)
#> Convert from m2 to km2
pcodes.clip$AREA_KM2_CLIP<- pcodes.clip$AreaM2 / 1000000
pcodes.clip$AREA_KM2_CLIP <- round(pcodes.clip$AREA_KM2_CLIP, digits = 3)
pcodes.clip$AREA_KM2_CLIP <- as.numeric(pcodes.clip$AREA_KM2_CLIP)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "pcodes.clip.shp"))
unlink(here("Out", "Tests", "pcodes.clip.dbf"))
unlink(here("Out", "Tests", "pcodes.clip.prj"))
unlink(here("Out", "Tests", "pcodes.clip.shx"))
st_write(pcodes.clip, here("Out", "Tests", "pcodes.clip.shp"), quiet = TRUE)


#> To calculate agri land differences for postcodes which straddle the boundary of target county...
#> 1. Intersect pcodes.int with corine.agri, to calculate area of agri land for all intersecting pcode district polys (whole polygon)
#> 2. Intersect pcodes.clip with corine.int to calculate area of agri land ony within target county polygon
#> 3. Calculate percentage difference between the TWO sets of area calculations to gove % of agri land in each pcode district polygon


#> 1. Intersect pcodes.int with corine.agri, to calculate area of agri land for all intersecting pcode district polys (whole polygon)
agri.int.total <- st_intersection(pcodes.int, corine.agri)
#> Calculate area of Corine agri land
agri.int.total$AREA_AGRI_TOTAL <- st_area(agri.int.total)
#> Convert from m2 to km2
agri.int.total$AREA_AGRI_TOTAL <- agri.int.total$AREA_AGRI_TOTAL / 1000000
agri.int.total$AREA_AGRI_TOTAL <- round(agri.int.total$AREA_AGRI_TOTAL, digits = 3)
agri.int.total$AREA_AGRI_TOTAL <- as.numeric(agri.int.total$AREA_AGRI_TOTAL)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "agri.int.total.shp"))
unlink(here("Out", "Tests", "agri.int.total.dbf"))
unlink(here("Out", "Tests", "agri.int.total.prj"))
unlink(here("Out", "Tests", "agri.int.total.shx"))
st_write(agri.int.total, here("Out", "Tests", "agri.int.total.shp"), quiet = TRUE)
#> Create grouped (by pcode) non-geo version of the table
st_geometry(agri.int.total) <- NULL
agri.int.total <- agri.int.total |> 
  select(name, AREA_AGRI_TOTAL) |> 
  group_by(name) |> 
  summarise(AGRI_TOTAL = sum(AREA_AGRI_TOTAL))


#> 2. Intersect pcodes.clip with corine.int to calculate area of agri land ony within target county polygon
agri.int.clip <- st_intersection(pcodes.clip, corine.agri)
#> Calculate area of Corine agri land
agri.int.clip$AREA_AGRI_CLIP <- st_area(agri.int.clip)
#> Convert from m2 to km2
agri.int.clip$AREA_AGRI_CLIP <- agri.int.clip$AREA_AGRI_CLIP / 1000000
agri.int.clip$AREA_AGRI_CLIP <- round(agri.int.clip$AREA_AGRI_CLIP, digits = 3)
agri.int.clip$AREA_AGRI_CLIP <- as.numeric(agri.int.clip$AREA_AGRI_CLIP)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "agri.int.clip.shp"))
unlink(here("Out", "Tests", "agri.int.clip.dbf"))
unlink(here("Out", "Tests", "agri.int.clip.prj"))
unlink(here("Out", "Tests", "agri.int.clip.shx"))
st_write(agri.int.clip, here("Out", "Tests", "agri.int.clip.shp"), quiet = TRUE)
#> Create grouped (by pcode) non-geo version of the table
st_geometry(agri.int.clip) <- NULL
agri.int.clip <- agri.int.clip |> 
  select(name, AREA_AGRI_CLIP) |> 
  group_by(name) |> 
  summarise(AGRI_CLIP = sum(AREA_AGRI_CLIP))


#> 3. Calculate percentage difference between the two sets of area calculations to gove % of agri land in each pcode district polygon
#> Merge the two data frames to show size of agri area within and outside target county
agri.merge <- merge(agri.int.total, agri.int.clip, by = "name", all.x = TRUE)
#> Round area calcs up to 1 decimal place (to negate effect of small amounts of agri land being lost when clipped)
agri.merge$AGRI_TOTAL <- round(agri.merge$AGRI_TOTAL, digits = 3)
agri.merge$AGRI_CLIP <- round(agri.merge$AGRI_CLIP, digits = 2)
#> Calculate percentage of agri land within target county 
agri.merge$PCENT_AGRI_LAND <- agri.merge$AGRI_CLIP/ agri.merge$AGRI_TOTAL
#> Round up
agri.merge$PCENT_AGRI_LAND <- round(agri.merge$PCENT_AGRI_LAND, digits = 3)



#> Extract postcode districts from the BPS source data ("db") which intersect with the active county polygon boundary
#> Do this via a merge between db and pcodes.int, keeping only matching records in a new data frame
#> Rename pcodes.in and drop geom
pcodes.active <- pcodes.int
st_geometry(pcodes.active) <- NULL
BPS.merge <- merge(db, pcodes.active, by.x = "PostcodePrefix_F202B", "name", all.x = FALSE)
#> Convert all NAs to zeros in data frame
BPS.merge[is.na(BPS.merge)] <- 0



#> BPS Calcs by year

#> 2020
#> Create table with baseline P1 (BPS + Greening) payments by postcode district for 2020 (no payments applied)
t1.2020 <- BPS.merge %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2020_Baseline = sum(`Basic payment scheme`) + sum(`Greening: practices beneficial for climate and environment`))
#> Check sum total of P1 for reference
sum.original.p1 <- sum(t1.2020$P1_2020_Baseline)
#> Merge with agri.merge data frame to append 2020 calcs and create new master table "t1"
t1 <- merge(t1.2020, agri.merge, by.x = "PostcodePrefix_F202B", by.y = "name")
#> #> Calculate P1 payment based on area of agri land within each postcode area
t1$P1_2020 <- t1$P1_2020_Baseline * t1$PCENT_AGRI_LAND
#> Round to 2 decimal places
t1$P1_2020 <- round(t1$P1_2020, digits = 2)



#> 2021
t1.2021 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2021$P1_2020 <- t1.2021$`Basic payment scheme` + t1.2021$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2021 payments
t1.2021 <- t1.2021 %>% 
  mutate(P1_2021 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.95,
                              (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.90,
                              (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.80,
                              P1_2020 > 150000 ~ P1_2020 * 0.75))
# Group by postcode an summarise
t1.2021 <- t1.2021 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2021 = sum(P1_2021))
#> Round
t1.2021$P1_2021 <- round(t1.2021$P1_2021, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2021, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2021 <- t1$P1_2021 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2021 <- round(t1$P1_2021, digits = 2)
#> Check sum
sum.p1.2021 <- sum(t1$P1_2021)



#> 2022
t1.2022 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2022$P1_2020 <- t1.2022$`Basic payment scheme` + t1.2022$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2022 payments
t1.2022 <- t1.2022 %>% 
   mutate(P1_2022 = case_when(P1_2020<=30000 ~ P1_2020* 0.80,
                             (P1_2020>30000 & P1_2020<=50000) ~ P1_2020* 0.75,
                             (P1_2020>50000 & P1_2020<=150000) ~ P1_2020* 0.65,
                             P1_2020> 150000 ~ P1_2020* 0.60))
# Group by postcode an summarise
t1.2022 <- t1.2022 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2022 = sum(P1_2022))
#> Round
t1.2022$P1_2022 <- round(t1.2022$P1_2022, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2022, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2022 <- t1$P1_2022 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2022 <- round(t1$P1_2022, digits = 2)
#> Check sum
sum.p1.2022 <- sum(t1$P1_2022)


#> 2023
t1.2023 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2023$P1_2020 <- t1.2023$`Basic payment scheme` + t1.2023$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2023 payments
t1.2023 <- t1.2023 %>% 
  mutate(P1_2023 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.65,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.60,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.50,
                             P1_2020 > 150000 ~ P1_2020 * 0.45))
# Group by postcode an summarise
t1.2023 <- t1.2023 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2023 = sum(P1_2023))
#> Round
t1.2023$P1_2023 <- round(t1.2023$P1_2023, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2023, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2023 <- t1$P1_2023 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2023 <- round(t1$P1_2023, digits = 2)
#> Check sum
sum.p1.2023 <- sum(t1$P1_2023)


#> 2024
t1.2024 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2024$P1_2020 <- t1.2024$`Basic payment scheme` + t1.2024$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2024 payments
t1.2024 <- t1.2024 %>% 
  mutate(P1_2024 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.50,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.45,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.35,
                             P1_2020 > 150000 ~ P1_2020 * 0.30))
# Group by postcode an summarise
t1.2024 <- t1.2024 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2024 = sum(P1_2024))
#> Round
t1.2024$P1_2024 <- round(t1.2024$P1_2024, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2024, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2024 <- t1$P1_2024 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2024 <- round(t1$P1_2024, digits = 2)
#> Check sum
sum.p1.2024 <- sum(t1$P1_2024)


#> 2025
t1.2025 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2025$P1_2020 <- t1.2025$`Basic payment scheme` + t1.2025$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2025 payments
t1.2025 <- t1.2025 %>% 
 mutate(P1_2025 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.40,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.35,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.25,
                             P1_2020 > 150000 ~ P1_2020 * 0.20))
# Group by postcode an summarise
t1.2025 <- t1.2025 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2025 = sum(P1_2025))
#> Round
t1.2025$P1_2025 <- round(t1.2025$P1_2025, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2025, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2025 <- t1$P1_2025 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2025 <- round(t1$P1_2025, digits = 2)
#> Check sum
sum.p1.2025 <- sum(t1$P1_2025)


#> 2026
t1.2026 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2026$P1_2020 <- t1.2026$`Basic payment scheme` + t1.2026$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2026 payments
t1.2026 <- t1.2026 %>% 
  mutate(P1_2026 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.25,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.25,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.15,
                             P1_2020 > 150000 ~ P1_2020 * 0.15))
# Group by postcode an summarise
t1.2026 <- t1.2026 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2026 = sum(P1_2026))
#> Round
t1.2026$P1_2026 <- round(t1.2026$P1_2026, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2026, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2026 <- t1$P1_2026 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2026 <- round(t1$P1_2026, digits = 2)
#> Check sum
sum.p1.2026 <- sum(t1$P1_2026)


#> 2027
t1.2027 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2027$P1_2020 <- t1.2027$`Basic payment scheme` + t1.2027$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2027 payments
t1.2027 <- t1.2027 %>% 
   mutate(P1_2027 = case_when(P1_2020<=30000 ~ P1_2020* 0.15,
                             (P1_2020>30000 & P1_2020<=50000) ~ P1_2020* 0.15,
                             (P1_2020>50000 & P1_2020<=150000) ~ P1_2020* 0.10,
                             P1_2020> 150000 ~ P1_2020* 0.10))
# Group by postcode an summarise
t1.2027 <- t1.2027 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2027 = sum(P1_2027))
#> Round
t1.2027$P1_2027 <- round(t1.2027$P1_2027, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2027, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2027 <- t1$P1_2027 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2027 <- round(t1$P1_2027, digits = 2)
#> Check sum
sum.p1.2027 <- sum(t1$P1_2027)

#> Convert all NAs to zeros in data frame
t1[is.na(t1)] <- 0



#> Create plot

#> P1 payments by year
p1.in <- as.data.frame(colSums(t1[6:13]))
#> Change rownames to column "Value"
p1.plot <- tibble::rownames_to_column(p1.in, "VALUE")
#> Extract year from string
p1.plot$VALUE <- sub("^.*([0-9]{4}).*", "\\1", p1.plot$VALUE)
names(p1.plot)[1]<-paste("Year")
names(p1.plot)[2]<-paste("Value_bps")

#> Print table t1 on web page
print(knitr::kable(t1, col.names = c("PC", "2020_Unadjusted", "Agri_Total", "Agri_In", "Agri_PC", "2020", "2021", "2022", "2023", "2024", "2025", "2026", "2027"), caption = paste0("Table 1.", i , " ", "P1 payments by postcode district - ",  active_county)))

#> Add line breaks
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")


# Simple bar chart
plot1 <- ggplot(p1.plot, aes(x=Year, y=Value_bps)) + 
  theme_bw() +
  geom_bar(stat = "identity", fill="#2c66b8", width = 0.6) +
  geom_text(aes(label= paste0(round(Value_bps / 1000000, digits = 1), " M")), vjust = 1.5, colour = "white") +
  labs(title = paste0("Pillar 1 payments:", " ", active_county), x = "Year", y = "Value (£)") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6))
plot1

#> Add line breaks
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")

#> Render on web page
print(plot1)

#> Export the plot
ggsave(here("Out", "P1", "County", "Plot", paste0("P1_payments_", active_county, ".png")))

#> Export final t1 table
write_csv(t1, here("Out", "P1", "County", "CSV", paste0("P1_payments_", active_county, ".csv")))

}
Table 1.1 P1 payments by postcode district - Cornwall
PC 2020_Unadjusted Agri_Total Agri_In Agri_PC 2020 2021 2022 2023 2024 2025 2026 2027
EX22 5669306.76 310.698 60.04 0.193 1094176.20 1004087.41 839960.98 675834.55 511708.12 402290.50 258047.81 156378.31
EX23 2746120.50 163.305 158.20 0.969 2660990.76 2454921.44 2055772.82 1656624.21 1257475.59 991376.52 637061.79 385055.67
EX39 4543093.48 221.849 0.38 0.002 9086.19 8161.25 6798.32 5435.39 4072.46 3163.85 2034.83 1244.57
PL10 189625.56 11.641 11.62 0.998 189246.31 161035.17 132648.23 104261.28 75874.34 56949.70 34812.37 22137.34
PL11 693638.88 32.302 32.29 1.000 693638.88 590828.92 486783.08 382737.25 278691.42 209327.53 140836.52 87759.23
PL12 1606968.13 97.248 96.39 0.991 1592505.42 1434640.46 1195764.66 956888.85 718013.04 558762.49 356392.72 218008.99
PL13 1346066.02 70.898 70.86 0.999 1344719.95 1223622.96 1021914.96 820206.97 618498.98 484026.99 317203.71 192219.86
PL14 3699213.10 225.157 225.16 1.000 3699213.10 3360619.05 2805737.08 2250855.12 1695973.15 1326051.84 855607.66 520284.16
PL15 6365797.99 366.707 332.65 0.907 5773778.78 5260089.03 4394022.21 3527955.40 2661888.58 2084510.70 1345868.83 817278.89
PL16 777462.94 51.306 0.40 0.008 6219.70 5753.88 4820.93 3887.97 2955.02 2333.05 1514.73 912.86
PL17 1330900.25 84.060 82.93 0.987 1313598.55 1205714.23 1008674.44 811634.66 614594.87 483235.02 311738.76 188709.35
PL18 149813.87 10.966 9.86 0.899 134682.67 127948.54 107746.14 87543.74 67341.34 53873.07 33670.67 20202.40
PL19 2872716.90 178.363 1.02 0.006 17236.30 15349.10 12763.66 10178.21 7592.77 5869.14 3750.89 2306.35
PL20 2201678.08 241.356 0.82 0.003 6605.03 5733.21 4742.45 3751.70 2760.94 2100.44 1344.66 837.45
PL22 1066944.05 51.812 51.81 1.000 1066944.05 975915.67 815874.06 655832.45 495790.85 389096.44 256291.14 154819.17
PL23 342779.53 19.564 19.53 0.998 342093.97 298884.07 247569.98 196255.88 144941.79 110732.39 71384.03 44244.36
PL24 295252.56 22.745 22.72 0.999 294957.31 278472.22 234228.63 189985.02 145741.43 116245.70 73739.33 44243.59
PL25 102286.28 9.102 9.10 1.000 102286.28 97171.97 81829.02 66486.08 51143.14 40914.51 25571.57 15342.94
PL26 2117478.59 116.673 116.64 1.000 2117478.59 1866270.68 1548648.89 1231027.11 913405.32 701657.46 441685.54 273779.74
PL27 2006347.61 105.911 105.88 1.000 2006347.61 1794473.10 1493520.96 1192568.82 891616.67 690981.91 442622.89 271470.14
PL28 395589.12 27.003 26.96 0.998 394797.94 357752.54 298532.86 239313.16 180093.47 140613.68 88187.79 53963.85
PL29 418500.38 17.702 17.67 0.998 417663.38 365148.36 302498.85 239849.34 177199.83 135433.50 84359.76 52621.47
PL30 4225344.65 246.323 246.32 1.000 4225344.65 3774675.54 3140873.85 2507072.15 1873270.45 1450735.99 921158.70 566212.97
PL31 123129.98 8.036 8.04 1.000 123129.98 114872.48 96402.98 77933.49 59463.99 47150.99 30782.50 18469.50
PL32 1462284.16 86.108 86.11 1.000 1462284.16 1289920.03 1070577.40 851234.78 631892.15 485663.74 309257.44 191185.82
PL33 367946.45 16.668 16.66 1.000 367946.45 341738.88 286546.91 231354.94 176162.98 139368.33 91986.61 55191.97
PL34 194194.17 14.574 14.54 0.998 193805.78 173623.95 144553.07 115482.21 86411.34 67030.76 41457.08 25573.68
PL35 405298.48 25.851 25.81 0.998 404487.88 378643.60 317970.41 257297.23 196624.05 156175.26 101121.97 60673.18
PL5 115480.60 6.203 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
TR1 150767.81 10.057 10.06 1.000 150767.81 141036.09 118420.92 95805.75 73190.58 58113.80 37691.95 22615.17
TR10 242695.35 23.471 23.47 1.000 242695.35 227005.59 190601.29 154196.98 117792.68 93523.15 60673.84 36404.30
TR11 809918.04 47.771 47.76 1.000 809918.04 739935.50 618447.79 496960.08 375472.38 294480.57 190239.64 115367.77
TR12 2347108.17 146.247 146.11 0.999 2344761.06 2114143.59 1762429.43 1410715.27 1059001.11 824525.00 527904.35 322571.20
TR13 1529143.44 78.919 78.90 1.000 1529143.44 1363864.31 1134492.79 905121.28 675749.76 522835.42 333307.45 204882.31
TR14 570914.31 52.522 52.52 1.000 570914.31 534880.93 449243.79 363606.64 277969.50 220878.06 142728.58 85637.15
TR15 52472.02 5.954 5.95 0.999 52419.55 49798.57 41935.64 34072.70 26209.77 20967.82 13104.89 7862.93
TR16 562427.41 59.317 59.31 1.000 562427.41 517473.13 433109.02 348744.91 264380.80 208138.06 129384.91 78753.14
TR17 31789.87 3.792 3.79 0.999 31758.08 30170.18 25406.47 20642.76 15879.04 12703.23 7939.52 4763.71
TR18 43280.41 2.995 3.00 1.002 43366.97 41198.62 34693.58 28188.53 21683.48 17346.78 10841.74 6505.04
TR19 1255959.76 84.766 84.40 0.996 1250935.92 1130289.08 942648.69 755008.30 567367.91 442274.32 281189.12 171867.96
TR2 3014939.68 158.572 158.50 1.000 3014939.68 2674622.68 2222381.73 1770140.77 1317899.82 1016405.85 650596.61 400671.80
TR20 1844232.47 107.401 107.38 1.000 1844232.47 1655324.13 1378689.26 1102054.39 825419.52 640996.27 401884.76 247048.19
TR21 112278.57 5.650 5.63 0.996 111829.46 97562.55 80788.14 64013.72 47239.30 36056.36 22173.75 13882.61
TR22 2863.23 1.720 1.69 0.983 2814.56 2673.83 2251.64 1829.46 1407.28 1125.82 703.64 422.18
TR23 0.00 1.689 1.68 0.995 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
TR24 15758.29 2.400 2.40 1.000 15758.29 14970.38 12606.63 10242.89 7879.15 6303.32 3939.57 2363.74
TR25 12436.24 2.839 2.78 0.979 12175.08 11566.33 9740.06 7913.81 6087.54 4870.04 3043.77 1826.27
TR26 509676.37 30.855 30.82 0.999 509166.69 462369.99 385994.99 309619.98 233244.98 182328.31 119584.25 72521.30
TR27 1128696.86 58.065 58.06 1.000 1128696.86 993178.89 823874.36 654569.83 485265.30 372395.61 235024.62 145729.73
TR3 777745.65 49.287 49.28 1.000 777745.65 703902.64 587240.79 470578.94 353917.10 276142.53 174426.84 106657.06
TR4 1695773.01 109.954 109.95 1.000 1695773.01 1533393.78 1279027.83 1024661.87 770295.92 600718.62 385095.57 234942.11
TR5 146924.33 15.864 15.86 1.000 146924.33 137698.29 115659.64 93620.99 71582.34 56889.91 36731.08 22038.65
TR6 114737.04 7.700 7.70 1.000 114737.04 101498.14 84287.59 67077.03 49866.48 38392.77 23682.90 14709.87
TR7 0.00 2.474 2.47 0.998 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
TR8 1918497.42 124.291 124.25 1.000 1918497.42 1702684.63 1414910.01 1127135.40 839360.79 647511.05 415615.20 255770.04
TR9 654565.36 42.333 42.33 1.000 654565.36 595145.15 496960.35 398775.54 300590.74 235134.20 149961.52 91344.89

Table 1.2 P1 payments by postcode district - Devon
PC 2020_Unadjusted Agri_Total Agri_In Agri_PC 2020 2021 2022 2023 2024 2025 2026 2027
DT6 3036994.48 158.966 0.36 0.002 6073.99 5382.63 4471.53 3560.43 2649.33 2041.93 1292.80 798.25
DT7 238968.44 16.753 15.32 0.914 218417.15 192242.36 159479.78 126717.21 93954.64 72112.92 45741.86 28331.36
EX1 194693.43 3.023 3.02 0.999 194498.74 183115.46 153940.65 124765.84 95591.03 76141.15 48624.69 29174.81
EX10 988503.32 55.885 55.88 1.000 988503.32 865341.81 717066.31 568790.81 420515.32 321664.98 211756.34 130590.75
EX11 515375.44 33.835 33.84 1.000 515375.44 464148.89 386842.58 309536.26 232229.94 180692.40 117277.18 71522.98
EX12 475378.22 24.567 24.56 1.000 475378.22 413573.15 342266.42 270959.68 199652.95 152115.13 97410.08 60589.49
EX13 2276377.07 111.911 107.43 0.960 2185321.99 2017742.73 1689944.43 1362146.13 1034347.83 815815.63 521702.15 315484.13
EX14 3584346.54 187.883 187.12 0.996 3570009.15 3288289.19 2752787.82 2217286.45 1681785.07 1324784.15 841566.24 510033.34
EX15 3044385.84 167.702 166.60 0.993 3023075.14 2826972.53 2373511.26 1920049.99 1466588.72 1164281.21 750243.74 450698.75
EX16 6005179.21 326.159 323.42 0.992 5957137.78 5548811.89 4655241.22 3761670.55 2868099.88 2272386.11 1443349.96 870603.42
EX17 5886022.85 293.242 293.24 1.000 5886022.85 5295893.63 4412990.20 3530086.78 2647183.35 2058581.06 1315300.89 804801.02
EX18 1746019.28 103.830 103.83 1.000 1746019.28 1618154.73 1356251.84 1094348.95 832446.06 657844.13 422280.55 254790.76
EX19 1737511.41 94.747 94.75 1.000 1737511.41 1624901.40 1364274.69 1103647.98 843021.27 669270.13 427002.09 256938.83
EX2 573837.75 19.150 19.15 1.000 573837.75 495976.30 409900.64 323824.98 237749.31 180365.54 112993.93 70842.91
EX20 7462904.40 443.972 443.97 1.000 7462904.40 6763333.19 5643897.53 4524461.87 3405026.21 2658735.77 1685950.37 1029547.80
EX21 2399858.70 146.828 146.83 1.000 2399858.70 2236945.26 1876966.45 1516987.65 1157008.84 917022.97 583169.41 351581.17
EX22 5669306.76 310.698 250.65 0.807 4575130.56 4198438.03 3512168.45 2825898.86 2139629.28 1682116.22 1078987.50 653872.02
EX23 2746120.50 163.305 5.11 0.031 85129.74 78537.22 65767.76 52998.30 40228.84 31715.86 20380.72 12318.60
EX24 781691.86 44.070 44.07 1.000 781691.86 728370.10 611116.32 493862.54 376608.76 298439.57 190347.70 114716.14
EX3 116221.08 7.142 7.14 1.000 116221.08 108799.55 91366.39 73933.23 56500.06 44877.96 29055.27 17433.16
EX31 4437568.50 217.860 217.86 1.000 4437568.50 4001220.85 3335585.57 2669950.30 2004315.02 1560558.17 1002774.60 612326.51
EX32 2503991.07 115.188 115.19 1.000 2503991.07 2258849.75 1883251.09 1507652.43 1132053.77 881654.66 562531.22 343865.39
EX33 721810.94 51.216 51.21 1.000 721810.94 643023.10 534751.46 426479.82 318208.18 246027.08 157286.94 96688.75
EX34 1182067.84 74.954 74.95 1.000 1182067.84 1094341.86 917031.69 739721.51 562411.34 444204.55 287578.34 173340.86
EX35 653966.01 64.790 54.03 0.834 545407.65 479120.50 397309.35 315498.20 233687.05 179146.29 112574.02 69922.20
EX36 4654512.01 230.731 220.83 0.957 4454367.99 4027568.80 3359413.60 2691258.40 2023103.21 1577666.41 1008182.35 615450.37
EX37 2107604.22 115.910 115.91 1.000 2107604.22 1961614.56 1645473.93 1329333.30 1013192.67 802432.24 519266.49 312323.35
EX38 1772073.96 94.594 94.59 1.000 1772073.96 1641940.91 1376129.82 1110318.72 844507.63 667300.23 426286.32 257445.01
EX39 4543093.48 221.849 221.29 0.997 4529464.20 4068382.06 3388962.43 2709542.79 2030123.16 1577176.74 1014360.67 620416.94
EX4 345514.09 33.458 33.46 1.000 345514.09 328238.39 276411.27 224584.16 172757.05 138205.64 86378.52 51827.11
EX5 3671497.74 208.646 208.65 1.000 3671497.74 3249242.29 2698517.62 2147792.96 1597068.30 1229918.53 776091.33 479833.11
EX6 2976600.96 173.859 173.86 1.000 2976600.96 2677064.48 2230574.33 1784084.19 1337594.04 1039933.95 674592.36 411711.21
EX7 246386.69 15.860 15.86 1.000 246386.69 234067.36 197109.35 160151.35 123193.34 98554.68 61596.67 36958.00
EX8 160334.56 13.778 13.78 1.000 160334.56 148533.85 124483.67 100433.48 76383.30 60349.84 40083.64 24050.18
EX9 735024.20 20.185 20.18 1.000 735024.20 599822.59 489568.96 379315.33 269061.70 195559.28 133165.17 84958.19
PL10 189625.56 11.641 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
PL12 1606968.13 97.248 0.86 0.009 14462.71 13029.03 10859.62 8690.21 6520.80 5074.53 3236.66 1979.90
PL15 6365797.99 366.707 34.06 0.093 592019.21 539347.61 450544.73 361741.84 272938.96 213737.04 137999.78 83800.37
PL16 777462.94 51.306 50.91 0.992 771243.24 713481.59 597795.10 482108.61 366422.12 289297.80 187827.07 113194.62
PL17 1330900.25 84.060 1.13 0.013 17301.70 15880.73 13285.48 10690.22 8094.97 6364.80 4105.98 2485.53
PL18 149813.87 10.966 1.10 0.100 14981.39 14232.32 11985.11 9737.90 7490.69 5992.56 3745.35 2247.21
PL19 2872716.90 178.363 177.34 0.994 2855480.60 2542834.90 2114512.81 1686190.72 1257868.63 972320.57 621397.02 382085.53
PL20 2201678.08 241.356 240.54 0.997 2195073.05 1905336.12 1576075.17 1246814.21 917553.26 698045.95 446874.13 278313.90
PL21 1668994.99 118.524 118.52 1.000 1668994.99 1528111.87 1277762.62 1027413.37 777064.12 610164.62 392921.22 238185.48
PL5 115480.60 6.203 6.20 1.000 115480.60 106244.85 88922.76 71600.67 54278.58 42730.52 28870.15 17322.09
PL6 358828.36 13.705 13.70 1.000 358828.36 297562.66 243738.41 189914.16 136089.90 100207.07 67211.32 42576.37
PL7 494378.82 47.727 47.73 1.000 494378.82 446623.20 372466.38 298309.56 224152.73 174714.85 110409.84 67564.39
PL8 1136984.10 50.334 50.22 0.998 1134710.13 968268.73 798062.22 627855.69 457649.18 344178.16 225296.70 141016.10
PL9 282611.66 19.426 19.36 0.997 281763.83 261252.70 218988.13 176723.55 134458.98 106282.60 70440.96 42264.57
TA20 2159808.09 117.094 2.78 0.024 51835.39 46462.20 38686.89 30911.58 23136.27 17952.73 11291.88 6941.82
TA21 1879680.39 96.095 22.65 0.236 443604.57 406890.73 340350.05 273809.36 207268.68 162908.22 104035.52 63107.87
TA22 1853621.93 96.166 16.65 0.173 320676.59 285914.84 237813.35 189711.86 141610.37 109542.71 70904.44 43469.13
TA24 3715605.73 308.095 11.38 0.037 137477.41 121563.89 100942.27 80320.66 59699.05 45951.31 29622.20 18248.04
TA3 2697932.22 168.419 2.53 0.015 40468.98 36520.47 30450.12 24379.78 18309.43 14262.53 9161.22 5592.34
TA4 4151861.21 262.718 2.65 0.010 41518.61 37824.49 31596.70 25368.91 19141.11 14989.25 9577.44 5826.69
TQ1 24333.43 2.753 2.75 0.999 24309.10 23093.64 19447.27 15800.91 12154.55 9723.64 6077.28 3646.36
TQ10 545741.00 57.013 57.01 1.000 545741.00 507933.39 426072.24 344211.09 262349.94 207775.84 131083.17 79185.11
TQ11 581560.02 49.024 49.02 1.000 581560.02 513472.29 426238.28 339004.28 251770.28 193614.28 123251.22 76164.61
TQ12 1454953.63 89.134 89.13 1.000 1454953.63 1333341.00 1115097.95 896854.91 678611.86 533116.50 344689.98 208718.83
TQ13 4300642.09 302.455 302.46 1.000 4300642.09 3904209.96 3259113.65 2614017.33 1968921.02 1538856.81 981892.72 598462.41
TQ14 153885.26 15.135 15.13 1.000 153885.26 144263.38 121180.59 98097.80 75015.01 59626.49 38471.32 23082.79
TQ2 24192.13 3.600 3.60 1.000 24192.13 22982.52 19353.70 15724.88 12096.07 9676.85 6048.03 3628.82
TQ3 216742.76 13.730 13.73 1.000 216742.76 192997.73 160486.32 127974.90 95463.49 73789.21 48255.77 29546.46
TQ4 97776.17 5.554 5.55 0.999 97678.39 89236.69 74584.93 59933.18 45281.41 35513.57 24419.60 14651.76
TQ5 139481.84 13.535 13.53 1.000 139481.84 122072.40 101150.13 80227.85 59305.57 45357.39 28992.11 17983.10
TQ6 688740.03 42.177 42.09 0.998 687362.55 610162.86 507058.48 403954.09 300849.71 232113.46 147258.76 90813.45
TQ7 2985679.51 176.474 176.18 0.998 2979708.15 2662548.61 2215592.39 1768636.16 1321679.94 1023709.13 661746.68 405366.04
TQ8 150169.89 10.807 10.74 0.994 149268.87 131730.69 109340.36 86950.03 64559.70 49632.82 32045.23 19754.34
TQ9 3454545.04 179.338 179.34 1.000 3454545.04 3106311.45 2588129.69 2069947.94 1551766.18 1206311.68 766868.52 469797.88

Table 1.3 P1 payments by postcode district - Dorset
PC 2020_Unadjusted Agri_Total Agri_In Agri_PC 2020 2021 2022 2023 2024 2025 2026 2027
BA12 6346310.03 295.497 0.66 0.002 12692.62 10276.15 8372.25 6468.36 4564.47 3295.21 2191.08 1412.86
BA21 237902.25 12.906 1.36 0.105 24979.74 23093.41 19346.45 15599.49 11852.53 9354.56 6244.93 3746.96
BA22 2812398.19 175.078 17.10 0.098 275615.02 247215.89 205873.64 164531.38 123189.13 95627.63 61442.81 37611.78
BA8 652094.50 38.458 1.62 0.042 27387.97 24637.10 20528.91 16420.71 12312.52 9573.72 6188.86 3779.13
BA9 1238362.79 67.125 0.68 0.010 12383.63 11033.29 9175.75 7318.20 5460.66 4222.30 2819.53 1719.36
BH10 18441.91 1.309 1.31 1.001 18460.35 17537.33 14768.28 11999.23 9230.18 7384.14 4615.09 2769.06
BH12 1679.04 0.161 0.16 0.994 1668.97 1585.52 1335.17 1084.83 834.48 667.59 417.24 250.35
BH15 1397.27 0.504 0.50 0.992 1386.09 1316.79 1108.88 900.96 693.05 554.44 346.53 207.91
BH16 214147.37 23.513 23.51 1.000 214147.37 199302.94 167180.84 135058.73 102936.63 81521.89 53536.84 32122.11
BH17 53650.86 3.205 3.20 0.998 53543.56 48530.75 40499.23 32467.69 24436.16 19081.80 13385.89 8031.53
BH18 24029.73 1.196 1.20 1.003 24101.82 22896.72 19281.45 15666.18 12050.90 9640.73 6025.45 3615.27
BH19 490262.12 35.884 35.84 0.999 489771.86 415393.98 341928.20 268462.42 194996.64 146019.45 90390.74 57439.66
BH20 3190961.34 176.773 176.71 1.000 3190961.34 2643555.09 2164910.89 1686266.69 1207622.49 888526.36 602366.94 380957.51
BH21 2516506.61 186.531 184.59 0.990 2491341.54 2072578.10 1698876.87 1325175.64 951474.41 702340.25 466047.66 295307.37
BH22 56651.95 9.775 9.78 1.001 56708.60 53873.17 45366.88 36860.59 28354.31 22683.44 14177.15 8506.29
BH23 600654.08 50.330 20.78 0.413 248070.14 208246.03 171035.51 133824.99 96614.47 71807.46 44264.99 28334.25
BH24 861956.35 90.496 12.01 0.133 114640.19 100938.70 83742.67 66546.64 49350.61 37886.59 24694.53 15213.27
BH25 250154.12 14.814 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
BH31 16081.54 4.238 3.94 0.930 14955.83 14208.04 11964.66 9721.29 7477.92 5982.34 3738.96 2243.37
BH6 7728.00 0.252 0.25 0.992 7666.18 7282.87 6132.94 4983.01 3833.09 3066.47 1916.54 1149.93
BH7 155716.17 0.168 0.17 1.012 157584.76 118188.58 94550.86 70913.15 47275.43 31516.95 23637.72 15758.48
BH8 40754.19 2.719 2.72 1.000 40754.19 36809.82 30696.69 24583.57 18470.44 14395.02 10188.55 6113.13
BH9 7150.73 1.055 1.05 0.995 7114.98 6759.22 5691.98 4624.73 3557.48 2845.99 1778.74 1067.25
DT1 321142.12 4.796 4.80 1.001 321463.26 260143.16 211923.67 163704.18 115484.69 83338.37 55307.14 35690.15
DT10 2215501.44 112.858 112.17 0.994 2202208.43 1965869.25 1635537.99 1305206.72 974875.46 754654.61 487547.94 298829.18
DT11 5268228.17 292.039 292.04 1.000 5268228.17 4387406.27 3597172.04 2806937.82 2016703.59 1489880.77 958786.31 611098.86
DT2 11058244.50 545.816 545.68 1.000 11058244.50 9217476.53 7558739.86 5900003.18 4241266.51 3135442.06 2049903.67 1301407.95
DT3 772281.56 83.856 83.86 1.000 772281.56 700608.06 584765.83 468923.59 353081.36 275853.20 178561.59 108587.83
DT4 98837.48 2.798 2.80 1.001 98936.32 82183.82 67343.38 52502.93 37662.47 27768.85 16863.63 10905.22
DT5 6850.52 4.455 4.44 0.997 6829.97 6488.47 5463.98 4439.48 3414.98 2731.99 1707.49 1024.50
DT6 3036994.48 158.966 158.59 0.998 3030920.49 2685929.95 2231291.88 1776653.80 1322015.73 1018923.68 645108.70 398327.36
DT7 238968.44 16.753 1.44 0.086 20551.29 18088.45 15005.76 11923.06 8840.37 6785.24 4303.94 2665.75
DT8 1105528.16 75.156 74.20 0.987 1091156.29 1001781.97 838108.53 674435.09 510761.64 401646.02 258597.52 156577.67
DT9 3850902.00 175.569 146.50 0.834 3211652.27 2788947.75 2307199.91 1825452.07 1343704.23 1022539.00 664634.44 412608.53
EX13 2276377.07 111.911 3.44 0.031 70567.69 65156.28 54571.12 43985.97 33400.82 26344.05 16846.63 10187.51
SP5 6376124.67 344.449 44.94 0.130 828896.21 678942.43 554608.00 430273.57 305939.14 223049.52 145190.22 93317.51
SP6 2077494.82 118.288 9.50 0.080 166199.59 139550.43 114620.49 89690.55 64760.62 48140.66 31784.20 20047.09
SP7 1987764.70 116.721 64.70 0.554 1101221.64 988172.61 822989.36 657806.12 492622.87 382500.71 246788.04 150924.56
SP8 1209961.34 73.093 70.14 0.960 1161562.89 1082177.87 907943.43 733709.00 559474.57 443318.28 284630.27 171354.20
TA18 1030727.57 51.620 4.59 0.089 91734.75 80285.49 66525.28 52765.07 39004.85 29831.38 19711.00 12148.87
TA20 2159808.09 117.094 19.12 0.163 352048.72 315555.75 262748.44 209941.13 157133.82 121928.95 76690.67 47146.55

Table 1.4 P1 payments by postcode district - Somerset
PC 2020_Unadjusted Agri_Total Agri_In Agri_PC 2020 2021 2022 2023 2024 2025 2026 2027
BA1 1317599.99 42.662 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
BA10 783558.26 43.613 43.51 0.998 781991.14 692668.34 575369.66 458070.99 340772.32 262573.21 165622.85 102361.20
BA11 2127620.07 115.923 114.02 0.984 2093578.15 1865741.39 1551704.67 1237667.95 923631.23 714273.41 455995.13 280337.02
BA12 6346310.03 295.497 1.15 0.004 25385.24 20552.30 16744.51 12936.72 9128.94 6590.41 4382.16 2825.71
BA13 820881.86 63.448 0.81 0.013 10671.46 9237.62 7636.90 6036.18 4435.46 3368.32 2145.41 1339.49
BA14 1248966.35 65.935 0.74 0.011 13738.63 12072.86 10012.06 7951.27 5890.47 4516.61 2843.26 1765.09
BA15 574707.38 31.390 0.09 0.003 1724.12 1486.77 1228.15 969.53 710.91 538.50 333.51 209.86
BA16 372739.20 23.243 23.24 1.000 372739.20 344620.29 288709.41 232798.53 176887.65 139613.73 88154.60 53395.78
BA2 1931256.06 120.317 19.45 0.162 312863.48 274609.56 227680.04 180750.51 133820.99 102534.64 65032.88 40338.03
BA20 64409.18 2.669 2.67 1.000 64409.18 61188.72 51527.34 41865.97 32204.59 25763.67 16102.30 9661.38
BA21 237902.25 12.906 11.55 0.895 212922.51 196843.86 164905.48 132967.10 101028.73 79736.47 53230.63 31938.38
BA22 2812398.19 175.078 157.98 0.902 2536783.17 2275395.23 1894877.75 1514360.29 1133842.81 880164.49 565524.61 346181.89
BA3 3760736.13 112.129 100.84 0.899 3380901.78 3022825.59 2515690.32 2008555.05 1501419.78 1163329.61 738957.05 454001.07
BA4 2751068.31 146.960 146.96 1.000 2751068.31 2415479.97 2002819.72 1590159.47 1177499.23 902392.40 586463.86 362008.64
BA5 1690939.49 117.359 117.36 1.000 1690939.49 1528378.10 1274737.17 1021096.25 767455.33 598361.38 382991.34 233769.16
BA6 1390647.45 83.723 83.72 1.000 1390647.45 1240386.27 1031789.15 823192.03 614594.91 475530.17 301497.93 185515.15
BA7 379069.44 24.908 24.91 1.000 379069.44 340820.01 283959.60 227099.18 170238.76 132331.82 87658.99 53306.23
BA8 652094.50 38.458 36.84 0.958 624706.53 561960.63 468254.65 374548.67 280842.69 218372.04 141165.04 86200.19
BA9 1238362.79 67.125 66.40 0.989 1224740.80 1091192.70 907481.58 723770.45 540059.34 417585.26 278851.66 170044.35
BS13 10404.83 2.755 1.69 0.613 6378.16 6059.25 5102.53 4145.80 3189.08 2551.26 1594.54 956.72
BS14 81151.54 6.419 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
BS15 10534.54 1.762 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
BS20 199815.22 17.663 17.66 1.000 199815.22 188009.81 158037.52 128065.24 98092.96 78111.44 49953.80 29972.28
BS21 548234.79 31.763 31.76 1.000 548234.79 505671.30 423436.08 341200.86 258965.64 204142.16 130747.15 79079.45
BS22 179628.40 14.406 14.41 1.000 179628.40 170646.98 143702.72 116758.46 89814.20 71851.36 44907.10 26944.26
BS23 2872.94 2.148 2.15 1.001 2875.81 2732.02 2300.65 1869.28 1437.91 1150.33 718.96 431.37
BS24 773092.66 46.867 46.87 1.000 773092.66 698733.71 582769.81 466805.91 350842.01 273532.74 173936.27 106295.45
BS25 262358.33 20.297 20.30 1.000 262358.33 247070.08 207716.33 168362.58 129008.83 102772.99 65589.58 39353.75
BS26 654617.58 44.695 44.70 1.000 654617.58 614658.27 516465.63 418273.00 320080.36 254618.60 163654.39 98192.64
BS27 315526.19 23.393 23.39 1.000 315526.19 298153.06 250824.14 203495.21 156166.28 124613.66 78881.55 47328.93
BS28 818646.31 46.399 46.40 1.000 818646.31 758089.58 635292.64 512495.69 389698.74 307834.11 198527.79 119730.05
BS29 136999.89 10.349 10.35 1.000 136999.89 130149.90 109599.91 89049.93 68499.94 54799.96 34249.97 20549.98
BS3 60737.92 1.009 0.47 0.466 28303.87 22730.02 18484.44 14238.86 9993.28 7162.89 4303.53 2859.36
BS30 318622.86 32.179 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
BS31 277742.52 17.951 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
BS39 804731.30 65.493 1.58 0.024 19313.55 17075.30 14178.26 11281.23 8384.20 6452.84 4219.54 2592.61
BS4 7065.85 0.694 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
BS40 1926969.72 126.762 80.36 0.634 1221698.80 1089277.65 906022.83 722768.01 539513.19 417343.31 269848.11 165466.53
BS41 302196.87 12.967 10.04 0.774 233900.38 210255.27 175170.22 140085.16 105000.10 81610.06 54588.47 33141.74
BS48 662720.47 39.377 39.38 1.000 662720.47 619520.94 520112.87 420704.80 321296.73 255024.68 165680.12 99408.07
BS49 445168.44 20.350 20.35 1.000 445168.44 414162.57 347387.30 280612.03 213836.77 169319.92 105460.48 63859.45
BS8 282824.38 11.121 11.12 1.000 282824.38 264712.79 222289.13 179865.47 137441.82 109159.38 70706.10 42423.66
BS9 14657.62 0.409 0.38 0.929 13616.93 12936.08 10893.55 8851.00 6808.46 5446.77 3404.23 2042.54
DT10 2215501.44 112.858 0.69 0.006 13293.01 11866.41 9872.46 7878.51 5884.56 4555.26 2942.95 1803.80
DT2 11058244.50 545.816 0.10 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
DT8 1105528.16 75.156 0.95 0.013 14371.87 13194.70 11038.92 8883.14 6727.36 5290.17 3406.05 2062.32
DT9 3850902.00 175.569 29.07 0.166 639249.73 555114.30 459226.84 363339.38 267451.92 203526.95 132289.35 82125.92
EX13 2276377.07 111.911 1.04 0.009 20487.39 18916.34 15843.23 12770.12 9697.01 7648.27 4890.96 2957.66
EX14 3584346.54 187.883 0.76 0.004 14337.39 13205.98 11055.37 8904.76 6754.16 5320.42 3379.78 2048.33
EX15 3044385.84 167.702 1.10 0.007 21310.70 19928.31 16731.70 13535.10 10338.49 8207.42 5288.73 3177.13
EX16 6005179.21 326.159 2.74 0.008 48041.43 44748.48 37542.27 30336.05 23129.84 18325.69 11639.92 7021.00
EX35 653966.01 64.790 10.76 0.166 108558.36 95364.51 79080.76 62797.00 46513.25 35657.41 22406.82 13917.37
EX36 4654512.01 230.731 9.90 0.043 200144.02 180967.04 150945.44 120923.84 90902.23 70887.83 45299.73 27653.47
SN13 864110.22 39.210 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
SN14 2795932.20 167.304 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
SP8 1209961.34 73.093 2.33 0.032 38718.76 36072.60 30264.78 24456.97 18649.15 14777.28 9487.68 5711.81
TA1 583543.77 3.439 3.44 1.000 583543.77 463194.84 375663.28 288131.71 200600.15 142245.77 97741.11 63459.15
TA10 2109692.92 99.423 99.42 1.000 2109692.92 1925927.87 1609473.93 1293019.99 976566.05 765596.76 495036.55 300260.60
TA11 1073232.78 75.066 75.07 1.000 1073232.78 967715.24 806730.32 645745.41 484760.49 377437.21 240518.34 147089.99
TA12 558660.27 27.257 27.26 1.000 558660.27 506390.68 422591.64 338792.60 254993.56 199127.54 127181.38 77557.20
TA13 460559.01 21.453 21.45 1.000 460559.01 412193.56 343109.71 274025.86 204942.00 158886.10 101043.32 62035.64
TA14 177904.61 10.714 10.71 1.000 177904.61 157532.63 130846.94 104161.25 77475.56 59685.10 38461.99 23678.61
TA15 29031.53 3.646 3.65 1.001 29060.56 27607.53 23248.45 18889.36 14530.29 11624.22 7265.14 4359.08
TA16 86162.14 4.651 4.65 1.000 86162.14 81854.03 68929.71 56005.39 43081.07 34464.86 21540.54 12924.32
TA17 178279.32 13.211 13.21 1.000 178279.32 147100.00 120358.11 93616.21 66874.31 49046.38 29726.26 19320.11
TA18 1030727.57 51.620 47.03 0.911 938992.82 821798.70 680949.78 540100.85 399251.93 305352.65 201760.92 124355.28
TA19 1411625.94 72.381 72.38 1.000 1411625.94 1240763.41 1029019.52 817275.63 605531.74 464369.14 301098.96 185840.13
TA2 818663.67 38.814 38.81 1.000 818663.67 731050.62 608251.07 485451.52 362651.97 280785.60 184573.91 112753.55
TA20 2159808.09 117.094 95.20 0.813 1755923.98 1573906.89 1310518.29 1047129.69 783741.10 608148.70 382512.35 235154.27
TA21 1879680.39 96.095 73.44 0.764 1436075.82 1317222.55 1101811.17 886399.81 670988.43 527380.85 336792.95 204298.37
TA22 1853621.93 96.166 79.51 0.827 1532945.34 1366772.09 1136830.29 906888.49 676946.69 523652.16 338947.80 207797.54
TA23 740420.68 48.768 48.77 1.000 740420.68 682778.47 571715.37 460652.27 349589.16 275547.10 177352.78 107186.91
TA24 3715605.73 308.095 296.71 0.963 3578128.32 3163946.53 2627227.29 2090508.04 1553788.79 1195975.95 770977.93 474942.17
TA3 2697932.22 168.419 165.89 0.985 2657463.24 2398177.56 1999558.08 1600938.59 1202319.10 936572.78 601587.10 367230.12
TA4 4151861.21 262.718 260.07 0.990 4110342.60 3744624.40 3128073.01 2511521.62 1894970.24 1483935.97 948166.86 576841.99
TA5 2925401.70 151.213 151.21 1.000 2925401.70 2579174.07 2140363.81 1701553.56 1262743.30 970203.13 634622.56 390446.32
TA6 733241.99 29.051 29.05 1.000 733241.99 671836.35 561850.06 451863.76 341877.46 268553.26 172014.72 104338.41
TA7 2616149.68 166.411 166.41 1.000 2616149.68 2393421.48 2000999.03 1608576.58 1216154.13 954539.16 617416.55 374112.02
TA8 229119.70 10.372 10.34 0.997 228432.34 208333.83 174068.97 139804.13 105539.27 82696.04 51323.49 31372.55
TA9 1067636.38 65.504 65.50 1.000 1067636.38 1007790.37 847644.91 687499.45 527354.00 420590.36 266909.10 160145.46

2. Pillar 1: District-level

Calculation of P1 payment payments at sub-county (district unitary authority) level. Uses same basic methodology as above, but the level of data management complexity is higher due to multiple area-within-area p1 reduction calculations.

#> Get district boundary data for SW counties
districts.sw <- st_read(here("In", "Shape", "Districts.SW.UPDATED.shp"), quiet = TRUE)
#> Convert Sw counties layer back to British National Grid
south.west <- st_transform(south.west, 27700)
#> Join districts with SW counties
districts.sw <- st_join(districts.sw, south.west, largest = TRUE) |> 
  filter(!is.na(NAME.y)) |> 
  select(NAME.x, NAME.y) |> 
  rename(NAME_DISTRICT = NAME.x, NAME_COUNTY = NAME.y) 

#? Remove string "(B)" from district names
# # districts.sw$NAME <- str_remove_all(districts.sw$NAME, "(B)")
districts.sw$NAME_DISTRICT <- gsub("\\s*\\([^\\)]+\\)", "", districts.sw$NAME_DISTRICT)
#> Trim WS
districts.sw$NAME_DISTRICT <- trimws(districts.sw$NAME_DISTRICT)
# head(districts.sw)

#> Create an empty data table to hold the P1 payments for districts - will be populated in below loop
district.summary <- data.frame(Name=character(),
                                P1_2020=as.numeric(),
                                P1_2021=as.numeric(),
                                P1_2022=as.numeric(), 
                                P1_2023=as.numeric(), 
                                P1_2024=as.numeric(), 
                                P1_2025=as.numeric(), 
                                P1_2026=as.numeric(), 
                                P1_2027=as.numeric())
                         


#> Sequential integer counter (for figure numbers)
i <-0

#> Initiate for loop
for(active_county in counties){

#> Loop counter (for table caption number)
i <- i+1


#> Initiate for loop here
#> For county in counties
target_county <- counties.sw |>
  filter(NAME == active_county)
head(target_county)


# #> Get the first county (example here without a loop = Cornwall)
# #> For county in counties
# target_county <- counties.gb |>
#   filter(NAME == "Cornwall")


#> Clip districts.sw to extent of counties (do for one county initially - Cornwall)
districts.clip <- st_intersection(districts.sw, target_county) |> 
#> Get only the districts that are officially within the active county
  #> (During the st_join process some slightly overlapping districts from other counties may be joined
  filter(NAME_COUNTY == active_county)

head(active_county)
head(districts.clip)

#> NEED TO START ANOTHER FOR LOOP HERE TO LOOP THROUGH DISTRICTS WITHIN COUNTY
#> FOR NOW CHOOSE ISLES OF SCILLY WITHOUT FOR LOOP
districts <- as.list(districts.clip$NAME_DISTRICT)

#> Initiate districts loop
for(district in districts){

active_district <- districts.clip |> 
  filter(NAME_DISTRICT == district)
head(active_district)

#> Select postcode district areas that intersect with active district polygon
pcodes.int <- pcodes.gb[active_district,]
head(pcodes.int)
#> Add area (km2) column
pcodes.int$AreaM2 <- st_area(pcodes.int)
#> Convert from m2 to km2
pcodes.int$AREA_KM2_TOTAL <- pcodes.int$AreaM2 / 1000000
pcodes.int$AREA_KM2_TOTAL <- round(pcodes.int$AREA_KM2_TOTAL, digits = 3)
pcodes.int$AREA_KM2_TOTAL <- as.numeric(pcodes.int$AREA_KM2_TOTAL)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "pcodes.int.shp"))
unlink(here("Out", "Tests", "pcodes.int.dbf"))
unlink(here("Out", "Tests", "pcodes.int.prj"))
unlink(here("Out", "Tests", "pcodes.int.shx"))
st_write(pcodes.int, here("Out", "Tests", "pcodes.int.shp"), quiet = TRUE)


#> Generate a layer of clipped (cookie cutter) postcodes using active district boundary
pcodes.clip <- st_intersection(pcodes.int, active_district)
#> Add area (km2) column
pcodes.clip$AreaM2 <- st_area(pcodes.clip)
#> Convert from m2 to km2
pcodes.clip$AREA_KM2_CLIP<- pcodes.clip$AreaM2 / 1000000
pcodes.clip$AREA_KM2_CLIP <- round(pcodes.clip$AREA_KM2_CLIP, digits = 3)
pcodes.clip$AREA_KM2_CLIP <- as.numeric(pcodes.clip$AREA_KM2_CLIP)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "pcodes.clip.shp"))
unlink(here("Out", "Tests", "pcodes.clip.dbf"))
unlink(here("Out", "Tests", "pcodes.clip.prj"))
unlink(here("Out", "Tests", "pcodes.clip.shx"))
st_write(pcodes.clip, here("Out", "Tests", "pcodes.clip.shp"), quiet = TRUE)


#> To calculate agri land differences for postcodes which straddle the boundary of district...
#> 1. Intersect pcodes.int with corine.agri, to calculate area of agri land for all intersecting pcode district polys (whole polygon)
#> 2. Intersect pcodes.clip with corine.int to calculate area of agri land ony within target district polygon
#> 3. Calculate percentage difference between the TWO sets of area calculations to gove % of agri land in each pcode district polygon

#> 1. Intersect pcodes.int with corine.agri, to calculate area of agri land for all intersecting pcode district polys (whole polygon)
agri.int.total <- st_intersection(pcodes.int, corine.agri)
head(agri.int.total)
#> Calculate area of Corine agri land
agri.int.total$AREA_AGRI_TOTAL <- st_area(agri.int.total)
#> Convert from m2 to km2
agri.int.total$AREA_AGRI_TOTAL <- agri.int.total$AREA_AGRI_TOTAL / 1000000
agri.int.total$AREA_AGRI_TOTAL <- round(agri.int.total$AREA_AGRI_TOTAL, digits = 3)
agri.int.total$AREA_AGRI_TOTAL <- as.numeric(agri.int.total$AREA_AGRI_TOTAL)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "agri.int.total.shp"))
unlink(here("Out", "Tests", "agri.int.total.dbf"))
unlink(here("Out", "Tests", "agri.int.total.prj"))
unlink(here("Out", "Tests", "agri.int.total.shx"))
st_write(agri.int.total, here("Out", "Tests", "agri.int.total.shp"), quiet = TRUE)
#> Create grouped (by pcode) non-geo version of the table
st_geometry(agri.int.total) <- NULL
agri.int.total <- agri.int.total |> 
  select(name, AREA_AGRI_TOTAL) |> 
  group_by(name) |> 
  summarise(AGRI_TOTAL = sum(AREA_AGRI_TOTAL))
head(agri.int.total)



#> 2. Intersect pcodes.clip with corine.int to calculate area of agri land ony within target district polygon
agri.int.clip <- st_intersection(pcodes.clip, corine.agri)
#> Calculate area of Corine agri land
agri.int.clip$AREA_AGRI_CLIP <- st_area(agri.int.clip)
#> Convert from m2 to km2
agri.int.clip$AREA_AGRI_CLIP <- agri.int.clip$AREA_AGRI_CLIP / 1000000
agri.int.clip$AREA_AGRI_CLIP <- round(agri.int.clip$AREA_AGRI_CLIP, digits = 3)
agri.int.clip$AREA_AGRI_CLIP <- as.numeric(agri.int.clip$AREA_AGRI_CLIP)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "agri.int.clip.shp"))
unlink(here("Out", "Tests", "agri.int.clip.dbf"))
unlink(here("Out", "Tests", "agri.int.clip.prj"))
unlink(here("Out", "Tests", "agri.int.clip.shx"))
st_write(agri.int.clip, here("Out", "Tests", "agri.int.clip.shp"), quiet = TRUE)
#> Create grouped (by pcode) non-geo version of the table
st_geometry(agri.int.clip) <- NULL
agri.int.clip <- agri.int.clip |> 
  select(name, AREA_AGRI_CLIP) |> 
  group_by(name) |> 
  summarise(AGRI_CLIP = sum(AREA_AGRI_CLIP))


#> 3. Calculate percentage difference between the two sets of area calculations to gove % of agri land in each pcode district polygon
#> Merge the two data frames to show size of agri area within and outside target county
agri.merge <- merge(agri.int.total, agri.int.clip, by = "name", all.x = TRUE)
#> Round area calcs up to 1 decimal place (to negate effect of small amounts of agri land being lost when clipped)
agri.merge$AGRI_TOTAL <- round(agri.merge$AGRI_TOTAL, digits = 3)
agri.merge$AGRI_CLIP <- round(agri.merge$AGRI_CLIP, digits = 2)
#> Calculate percentage of agri land within target county 
agri.merge$PCENT_AGRI_LAND <- agri.merge$AGRI_CLIP/ agri.merge$AGRI_TOTAL
#> Round up
agri.merge$PCENT_AGRI_LAND <- round(agri.merge$PCENT_AGRI_LAND, digits = 3)
head(agri.merge)


#> Extract postcode districts from the BPS source data ("db") which intersect with the active county polygon boundary
#> Do this via a merge between db and pcodes.int, keeping only matching records in a new data frame
#> Rename pcodes.in and drop geom
pcodes.active <- pcodes.int
st_geometry(pcodes.active) <- NULL
BPS.merge <- merge(db, pcodes.active, by.x = "PostcodePrefix_F202B", "name", all.x = FALSE)
#> Convert all NAs to zeros in data frame
BPS.merge[is.na(BPS.merge)] <- 0
head(BPS.merge)


#> BPS Calcs by year (for target district)

#> 2020
#> Create table with baseline P1 (BPS + Greening) payments by postcode district for 2020 (no payments applied)
t1.2020 <- BPS.merge %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2020_Baseline = sum(`Basic payment scheme`) + sum(`Greening: practices beneficial for climate and environment`))
#> Check sum total of P1 for reference
sum.original.p1 <- sum(t1.2020$P1_2020_Baseline)
#> Merge with agri.merge data frame to append 2020 calcs and create new master table "t1"
t1 <- merge(t1.2020, agri.merge, by.x = "PostcodePrefix_F202B", by.y = "name")
#> #> Calculate P1 payment based on area of agri land within each postcode area
t1$P1_2020 <- t1$P1_2020_Baseline * t1$PCENT_AGRI_LAND
#> Round to 2 decimal places
t1$P1_2020 <- round(t1$P1_2020, digits = 2)




#> 2021
t1.2021 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2021$P1_2020 <- t1.2021$`Basic payment scheme` + t1.2021$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2021 payments
t1.2021 <- t1.2021 %>% 
  mutate(P1_2021 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.95,
                              (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.90,
                              (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.80,
                              P1_2020 > 150000 ~ P1_2020 * 0.75))
# Group by postcode an summarise
t1.2021 <- t1.2021 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2021 = sum(P1_2021))
#> Round
t1.2021$P1_2021 <- round(t1.2021$P1_2021, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2021, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2021 <- t1$P1_2021 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2021 <- round(t1$P1_2021, digits = 2)
#> Check sum
sum.p1.2021 <- sum(t1$P1_2021)



#> 2022
t1.2022 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2022$P1_2020 <- t1.2022$`Basic payment scheme` + t1.2022$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2022 payments
t1.2022 <- t1.2022 %>% 
   mutate(P1_2022 = case_when(P1_2020<=30000 ~ P1_2020* 0.80,
                             (P1_2020>30000 & P1_2020<=50000) ~ P1_2020* 0.75,
                             (P1_2020>50000 & P1_2020<=150000) ~ P1_2020* 0.65,
                             P1_2020> 150000 ~ P1_2020* 0.60))
# Group by postcode an summarise
t1.2022 <- t1.2022 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2022 = sum(P1_2022))
#> Round
t1.2022$P1_2022 <- round(t1.2022$P1_2022, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2022, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2022 <- t1$P1_2022 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2022 <- round(t1$P1_2022, digits = 2)
#> Check sum
sum.p1.2022 <- sum(t1$P1_2022)


#> 2023
t1.2023 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2023$P1_2020 <- t1.2023$`Basic payment scheme` + t1.2023$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2023 payments
t1.2023 <- t1.2023 %>% 
  mutate(P1_2023 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.65,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.60,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.50,
                             P1_2020 > 150000 ~ P1_2020 * 0.45))
# Group by postcode an summarise
t1.2023 <- t1.2023 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2023 = sum(P1_2023))
#> Round
t1.2023$P1_2023 <- round(t1.2023$P1_2023, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2023, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2023 <- t1$P1_2023 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2023 <- round(t1$P1_2023, digits = 2)
#> Check sum
sum.p1.2023 <- sum(t1$P1_2023)


#> 2024
t1.2024 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2024$P1_2020 <- t1.2024$`Basic payment scheme` + t1.2024$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2024 payments
t1.2024 <- t1.2024 %>% 
  mutate(P1_2024 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.50,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.45,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.35,
                             P1_2020 > 150000 ~ P1_2020 * 0.30))
# Group by postcode an summarise
t1.2024 <- t1.2024 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2024 = sum(P1_2024))
#> Round
t1.2024$P1_2024 <- round(t1.2024$P1_2024, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2024, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2024 <- t1$P1_2024 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2024 <- round(t1$P1_2024, digits = 2)
#> Check sum
sum.p1.2024 <- sum(t1$P1_2024)


#> 2025
t1.2025 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2025$P1_2020 <- t1.2025$`Basic payment scheme` + t1.2025$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2025 payments
t1.2025 <- t1.2025 %>% 
 mutate(P1_2025 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.40,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.35,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.25,
                             P1_2020 > 150000 ~ P1_2020 * 0.20))
# Group by postcode an summarise
t1.2025 <- t1.2025 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2025 = sum(P1_2025))
#> Round
t1.2025$P1_2025 <- round(t1.2025$P1_2025, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2025, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2025 <- t1$P1_2025 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2025 <- round(t1$P1_2025, digits = 2)
#> Check sum
sum.p1.2025 <- sum(t1$P1_2025)


#> 2026
t1.2026 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2026$P1_2020 <- t1.2026$`Basic payment scheme` + t1.2026$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2026 payments
t1.2026 <- t1.2026 %>% 
  mutate(P1_2026 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.25,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.25,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.15,
                             P1_2020 > 150000 ~ P1_2020 * 0.15))
# Group by postcode an summarise
t1.2026 <- t1.2026 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2026 = sum(P1_2026))
#> Round
t1.2026$P1_2026 <- round(t1.2026$P1_2026, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2026, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2026 <- t1$P1_2026 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2026 <- round(t1$P1_2026, digits = 2)
#> Check sum
sum.p1.2026 <- sum(t1$P1_2026)


#> 2027
t1.2027 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2027$P1_2020 <- t1.2027$`Basic payment scheme` + t1.2027$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2027 payments
t1.2027 <- t1.2027 %>% 
   mutate(P1_2027 = case_when(P1_2020<=30000 ~ P1_2020* 0.15,
                             (P1_2020>30000 & P1_2020<=50000) ~ P1_2020* 0.15,
                             (P1_2020>50000 & P1_2020<=150000) ~ P1_2020* 0.10,
                             P1_2020> 150000 ~ P1_2020* 0.10))
# Group by postcode an summarise
t1.2027 <- t1.2027 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2027 = sum(P1_2027))
#> Round
t1.2027$P1_2027 <- round(t1.2027$P1_2027, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2027, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2027 <- t1$P1_2027 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2027 <- round(t1$P1_2027, digits = 2)
#> Check sum
sum.p1.2027 <- sum(t1$P1_2027)

#> Convert all NAs to zeros in data frame
t1[is.na(t1)] <- 0

#> Create a new data frame to hold P1 reduction values by year summed for each district
head(active_district)
district.summary <- district.summary |> 
  add_row(Name = district, 
          P1_2020 = sum(t1$P1_2020),
          P1_2021 = sum(t1$P1_2021),
          P1_2022 = sum(t1$P1_2022),
          P1_2023 = sum(t1$P1_2023),
          P1_2024 = sum(t1$P1_2024),
          P1_2025 = sum(t1$P1_2025),
          P1_2026 = sum(t1$P1_2026),
          P1_2027 = sum(t1$P1_2027))
head(district.summary)
#> Wipe t1 rows clean
t1 <- t1[0,]



#> Web plot

}



#> Print table district.summary on web page
print(knitr::kable(district.summary, caption = paste0("Table 2.", i , " ", "P1 payments by district - ",  active_county)))

cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")

#> Create bar plots

#> Get a list of districts from the active district.summary table
districts.list <- as.list(district.summary$Name)

for (district.active in districts.list)

{

#> Create data table with only active row (district) from districts.summary table
df <- district.summary|> 
  filter(Name == district.active)
#> Shape data
p2.in <- as.data.frame(colSums(df[2:8]))
#> Change rownames to column "Value"
p2.plot <- tibble::rownames_to_column(p2.in, "VALUE")
#> Extract year from string
p2.plot$VALUE <- sub("^.*([0-9]{4}).*", "\\1", p2.plot$VALUE)
names(p2.plot)[1]<-paste("Year")
names(p2.plot)[2]<-paste("Value_bps")

# Simple bar chart
plot2 <- ggplot(p2.plot, aes(x=Year, y=Value_bps)) + 
  theme_bw() +
  geom_bar(stat = "identity", fill="#9e66ab", width = 0.6) +
  geom_text(aes(label= paste0(round(Value_bps / 1000000, digits = 2), " M")), vjust = 1.5, colour = "white") +
  labs(title = paste0("Pillar 1 payments:", " ", district.active), x = "Year", y = "Value (£)") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6))
plot2
print(plot2)


#> Export the plot
ggsave(here("Out", "P1", "District", "Plot", paste0("P1_payments_", district.active, ".png")))


}



#> Export table
write_csv(district.summary, here("Out", "P1", "District", "CSV", paste0("P1_payments_", active_county, ".csv")))



#> Wipe district.summary rows clean
district.summary <- district.summary[0,]


  
}
Table 2.1 P1 payments by district - Cornwall
Name P1_2020 P1_2021 P1_2022 P1_2023 P1_2024 P1_2025 P1_2026 P1_2027
Cornwall 51444463.2 46411032.5 38694363.0 30977693.47 23261024.03 18116577.69 11605951.87 7089087.51
Isles of Scilly 142447.9 126658.8 105291.6 83924.44 62557.25 48312.46 29834.18 18478.28

Table 2.2 P1 payments by district - Devon
Name P1_2020 P1_2021 P1_2022 P1_2023 P1_2024 P1_2025 P1_2026 P1_2027
Torbay 259792.3 234261.77 195292.91 156324.06 117355.23 91375.99 60112.04 36550.84
City of Plymouth 100824.6 91521.45 76397.76 61274.06 46150.36 36067.90 23602.91 14322.06
North Devon District 18008097.5 16347562.25 13646347.61 10945132.98 8243918.36 6443108.60 4129816.42 2515110.66
East Devon District 12996715.1 11742429.24 9792921.96 7843414.68 5893907.40 4594235.88 2942592.18 1796213.93
Teignbridge District 7934517.8 7191798.84 6001621.15 4811443.51 3621265.82 2827814.06 1812910.51 1104818.20
West Devon District 15685530.4 14127177.23 11774347.68 9421518.12 7068688.56 5500135.53 3505846.54 2145061.54
Mid Devon District 15601533.4 14292197.56 11951967.55 9611737.55 7271507.53 5711354.18 3646422.34 2213249.50
Exeter District 187590.0 172869.29 144730.80 116592.30 88453.80 69694.80 43963.44 26671.46
South Hams District 12558179.7 11239050.35 9355323.38 7471596.42 5587869.47 4332051.51 2783520.00 1705714.48
Torridge District 15909701.2 14575552.02 12189096.86 9802641.67 7416186.50 5825216.36 3733048.91 2264266.98

Table 2.3 P1 payments by district - Dorset
Name P1_2020 P1_2021 P1_2022 P1_2023 P1_2024 P1_2025 P1_2026 P1_2027
Bournemouth, Christchurch and Poole 620515.8 517037.9 423960.5 330883.2 237805.8 175754.2 116739.1 73882.43
Dorset 37717075.0 32230428.8 26572867.5 20915306.2 15257745.0 11486037.5 7466069.0 4675961.40

Table 2.4 P1 payments by district - Somerset
Name P1_2020 P1_2021 P1_2022 P1_2023 P1_2024 P1_2025 P1_2026 P1_2027
North Somerset 4326730 3981035 3332026 2683016 2034007 1601334 1026217 621277
Mendip District 12844498 11459629 9532954 7606280 5679605 4395155 2811340 1726783
Somerset West and Taunton District 16407836 14706999 12245824 9784648 7323473 5682689 3661983 2241188
South Somerset District 15587826 13954474 11616300 9278126 6939952 5381169 3469378 2124385
Sedgemoor District 8435821 7705442 6440069 5174695 3909322 3065740 1974880 1198335

3. Pillar 1: National Parks

Pillar one reduction analysis conducted for Exmoor and Dartmoor National Parks. National Parks boundary data from Natural England Geoportal.

#> Import nat parks data
nat.parks <- st_read(here("In", "Shape", "Nat.Parks.SW.shp"), stringsAsFactors = FALSE, quiet = TRUE)
#> Change nat parks "NAME" to Capital case
nat.parks <- nat.parks |> 
  mutate(NAME = str_to_title(NAME, locale = "en"))
# glimpse(nat.parks)

#> Create filter string (get names first!)
nat.parks.list <- c("Exmoor", "Dartmoor")

#> Sequential integer counter (for figure numbers)
i <-0 

#> Initiate loop
for(active_park in nat.parks.list){
  
#> Loop counter (for table caption number)
i <- i+1

#> Select target county (this will be start of loop)
park <- nat.parks|> 
  filter(NAME == active_park)

#> Select postcode district areas that intersect with active county polygon
pcodes.int <- pcodes.gb[park,]
#> Add area (km2) column
pcodes.int$AreaM2 <- st_area(pcodes.int)
#> Convert from m2 to km2
pcodes.int$AREA_KM2_TOTAL <- pcodes.int$AreaM2 / 1000000
pcodes.int$AREA_KM2_TOTAL <- round(pcodes.int$AREA_KM2_TOTAL, digits = 3)
pcodes.int$AREA_KM2_TOTAL <- as.numeric(pcodes.int$AREA_KM2_TOTAL)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "pcodes.int.shp"))
unlink(here("Out", "Tests", "pcodes.int.dbf"))
unlink(here("Out", "Tests", "pcodes.int.prj"))
unlink(here("Out", "Tests", "pcodes.int.shx"))
st_write(pcodes.int, here("Out", "Tests", "pcodes.int.shp"), quiet = TRUE)


#> Generate a layer of clipped (cookie cutter) postcodes using active county boundary
pcodes.clip <- st_intersection(pcodes.int, park)
#> Add area (km2) column
pcodes.clip$AreaM2 <- st_area(pcodes.clip)
#> Convert from m2 to km2
pcodes.clip$AREA_KM2_CLIP<- pcodes.clip$AreaM2 / 1000000
pcodes.clip$AREA_KM2_CLIP <- round(pcodes.clip$AREA_KM2_CLIP, digits = 3)
pcodes.clip$AREA_KM2_CLIP <- as.numeric(pcodes.clip$AREA_KM2_CLIP)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "pcodes.clip.shp"))
unlink(here("Out", "Tests", "pcodes.clip.dbf"))
unlink(here("Out", "Tests", "pcodes.clip.prj"))
unlink(here("Out", "Tests", "pcodes.clip.shx"))
st_write(pcodes.clip, here("Out", "Tests", "pcodes.clip.shp"), quiet = TRUE)


#> To calculate agri land differences for postcodes which straddle the boundary of target county...
#> 1. Intersect pcodes.int with corine.agri, to calculate area of agri land for all intersecting pcode district polys (whole polygon)
#> 2. Intersect pcodes.clip with corine.int to calculate area of agri land ony within target county polygon
#> 3. Calculate percentage difference between the TWO sets of area calculations to gove % of agri land in each pcode district polygon


#> 1. Intersect pcodes.int with corine.agri, to calculate area of agri land for all intersecting pcode district polys (whole polygon)
agri.int.total <- st_intersection(pcodes.int, corine.agri)
#> Calculate area of Corine agri land
agri.int.total$AREA_AGRI_TOTAL <- st_area(agri.int.total)
#> Convert from m2 to km2
agri.int.total$AREA_AGRI_TOTAL <- agri.int.total$AREA_AGRI_TOTAL / 1000000
agri.int.total$AREA_AGRI_TOTAL <- round(agri.int.total$AREA_AGRI_TOTAL, digits = 3)
agri.int.total$AREA_AGRI_TOTAL <- as.numeric(agri.int.total$AREA_AGRI_TOTAL)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "agri.int.total.shp"))
unlink(here("Out", "Tests", "agri.int.total.dbf"))
unlink(here("Out", "Tests", "agri.int.total.prj"))
unlink(here("Out", "Tests", "agri.int.total.shx"))
st_write(agri.int.total, here("Out", "Tests", "agri.int.total.shp"), quiet = TRUE)
#> Create grouped (by pcode) non-geo version of the table
st_geometry(agri.int.total) <- NULL
agri.int.total <- agri.int.total |> 
  select(name, AREA_AGRI_TOTAL) |> 
  group_by(name) |> 
  summarise(AGRI_TOTAL = sum(AREA_AGRI_TOTAL))


#> 2. Intersect pcodes.clip with corine.int to calculate area of agri land ony within target county polygon
agri.int.clip <- st_intersection(pcodes.clip, corine.agri)
#> Calculate area of Corine agri land
agri.int.clip$AREA_AGRI_CLIP <- st_area(agri.int.clip)
#> Convert from m2 to km2
agri.int.clip$AREA_AGRI_CLIP <- agri.int.clip$AREA_AGRI_CLIP / 1000000
agri.int.clip$AREA_AGRI_CLIP <- round(agri.int.clip$AREA_AGRI_CLIP, digits = 3)
agri.int.clip$AREA_AGRI_CLIP <- as.numeric(agri.int.clip$AREA_AGRI_CLIP)
#> Export data for detailed visual checking in QGIS
#> Remove existing shapefile
unlink(here("Out", "Tests", "agri.int.clip.shp"))
unlink(here("Out", "Tests", "agri.int.clip.dbf"))
unlink(here("Out", "Tests", "agri.int.clip.prj"))
unlink(here("Out", "Tests", "agri.int.clip.shx"))
st_write(agri.int.clip, here("Out", "Tests", "agri.int.clip.shp"), quiet = TRUE)
#> Create grouped (by pcode) non-geo version of the table
st_geometry(agri.int.clip) <- NULL
agri.int.clip <- agri.int.clip |> 
  select(name, AREA_AGRI_CLIP) |> 
  group_by(name) |> 
  summarise(AGRI_CLIP = sum(AREA_AGRI_CLIP))


#> 3. Calculate percentage difference between the two sets of area calculations to gove % of agri land in each pcode district polygon
#> Merge the two data frames to show size of agri area within and outside target county
agri.merge <- merge(agri.int.total, agri.int.clip, by = "name", all.x = TRUE)
#> Round area calcs up to 1 decimal place (to negate effect of small amounts of agri land being lost when clipped)
agri.merge$AGRI_TOTAL <- round(agri.merge$AGRI_TOTAL, digits = 3)
agri.merge$AGRI_CLIP <- round(agri.merge$AGRI_CLIP, digits = 2)
#> Calculate percentage of agri land within target county 
agri.merge$PCENT_AGRI_LAND <- agri.merge$AGRI_CLIP/ agri.merge$AGRI_TOTAL
#> Round up
agri.merge$PCENT_AGRI_LAND <- round(agri.merge$PCENT_AGRI_LAND, digits = 3)



#> Extract postcode districts from the BPS source data ("db") which intersect with the active county polygon boundary
#> Do this via a merge between db and pcodes.int, keeping only matching records in a new data frame
#> Rename pcodes.in and drop geom
pcodes.active <- pcodes.int
st_geometry(pcodes.active) <- NULL
BPS.merge <- merge(db, pcodes.active, by.x = "PostcodePrefix_F202B", "name", all.x = FALSE)
#> Convert all NAs to zeros in data frame
BPS.merge[is.na(BPS.merge)] <- 0



#> BPS Calcs by year

#> 2020
#> Create table with baseline P1 (BPS + Greening) payments by postcode district for 2020 (no payments applied)
t1.2020 <- BPS.merge %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2020_Baseline = sum(`Basic payment scheme`) + sum(`Greening: practices beneficial for climate and environment`))
#> Check sum total of P1 for reference
sum.original.p1 <- sum(t1.2020$P1_2020_Baseline)
#> Merge with agri.merge data frame to append 2020 calcs and create new master table "t1"
t1 <- merge(t1.2020, agri.merge, by.x = "PostcodePrefix_F202B", by.y = "name")
#> #> Calculate P1 payment based on area of agri land within each postcode area
t1$P1_2020 <- t1$P1_2020_Baseline * t1$PCENT_AGRI_LAND
#> Round to 2 decimal places
t1$P1_2020 <- round(t1$P1_2020, digits = 2)



#> 2021
t1.2021 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2021$P1_2020 <- t1.2021$`Basic payment scheme` + t1.2021$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2021 payments
t1.2021 <- t1.2021 %>% 
  mutate(P1_2021 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.95,
                              (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.90,
                              (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.80,
                              P1_2020 > 150000 ~ P1_2020 * 0.75))
# Group by postcode an summarise
t1.2021 <- t1.2021 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2021 = sum(P1_2021))
#> Round
t1.2021$P1_2021 <- round(t1.2021$P1_2021, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2021, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2021 <- t1$P1_2021 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2021 <- round(t1$P1_2021, digits = 2)
#> Check sum
sum.p1.2021 <- sum(t1$P1_2021)



#> 2022
t1.2022 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2022$P1_2020 <- t1.2022$`Basic payment scheme` + t1.2022$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2022 payments
t1.2022 <- t1.2022 %>% 
   mutate(P1_2022 = case_when(P1_2020<=30000 ~ P1_2020* 0.80,
                             (P1_2020>30000 & P1_2020<=50000) ~ P1_2020* 0.75,
                             (P1_2020>50000 & P1_2020<=150000) ~ P1_2020* 0.65,
                             P1_2020> 150000 ~ P1_2020* 0.60))
# Group by postcode an summarise
t1.2022 <- t1.2022 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2022 = sum(P1_2022))
#> Round
t1.2022$P1_2022 <- round(t1.2022$P1_2022, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2022, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2022 <- t1$P1_2022 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2022 <- round(t1$P1_2022, digits = 2)
#> Check sum
sum.p1.2022 <- sum(t1$P1_2022)


#> 2023
t1.2023 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2023$P1_2020 <- t1.2023$`Basic payment scheme` + t1.2023$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2023 payments
t1.2023 <- t1.2023 %>% 
  mutate(P1_2023 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.65,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.60,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.50,
                             P1_2020 > 150000 ~ P1_2020 * 0.45))
# Group by postcode an summarise
t1.2023 <- t1.2023 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2023 = sum(P1_2023))
#> Round
t1.2023$P1_2023 <- round(t1.2023$P1_2023, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2023, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2023 <- t1$P1_2023 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2023 <- round(t1$P1_2023, digits = 2)
#> Check sum
sum.p1.2023 <- sum(t1$P1_2023)


#> 2024
t1.2024 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2024$P1_2020 <- t1.2024$`Basic payment scheme` + t1.2024$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2024 payments
t1.2024 <- t1.2024 %>% 
  mutate(P1_2024 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.50,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.45,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.35,
                             P1_2020 > 150000 ~ P1_2020 * 0.30))
# Group by postcode an summarise
t1.2024 <- t1.2024 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2024 = sum(P1_2024))
#> Round
t1.2024$P1_2024 <- round(t1.2024$P1_2024, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2024, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2024 <- t1$P1_2024 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2024 <- round(t1$P1_2024, digits = 2)
#> Check sum
sum.p1.2024 <- sum(t1$P1_2024)


#> 2025
t1.2025 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2025$P1_2020 <- t1.2025$`Basic payment scheme` + t1.2025$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2025 payments
t1.2025 <- t1.2025 %>% 
 mutate(P1_2025 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.40,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.35,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.25,
                             P1_2020 > 150000 ~ P1_2020 * 0.20))
# Group by postcode an summarise
t1.2025 <- t1.2025 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2025 = sum(P1_2025))
#> Round
t1.2025$P1_2025 <- round(t1.2025$P1_2025, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2025, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2025 <- t1$P1_2025 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2025 <- round(t1$P1_2025, digits = 2)
#> Check sum
sum.p1.2025 <- sum(t1$P1_2025)


#> 2026
t1.2026 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2026$P1_2020 <- t1.2026$`Basic payment scheme` + t1.2026$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2026 payments
t1.2026 <- t1.2026 %>% 
  mutate(P1_2026 = case_when(P1_2020 <=30000 ~ P1_2020 * 0.25,
                             (P1_2020 >30000 & P1_2020 <=50000) ~ P1_2020 * 0.25,
                             (P1_2020 >50000 & P1_2020 <=150000) ~ P1_2020 * 0.15,
                             P1_2020 > 150000 ~ P1_2020 * 0.15))
# Group by postcode an summarise
t1.2026 <- t1.2026 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2026 = sum(P1_2026))
#> Round
t1.2026$P1_2026 <- round(t1.2026$P1_2026, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2026, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2026 <- t1$P1_2026 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2026 <- round(t1$P1_2026, digits = 2)
#> Check sum
sum.p1.2026 <- sum(t1$P1_2026)


#> 2027
t1.2027 <- BPS.merge %>% 
  select(PostcodePrefix_F202B, `Basic payment scheme`, `Greening: practices beneficial for climate and environment`)
#> Column to hold total total P1 payments (i.e. BPS + Greening) - baseline for 2020
t1.2027$P1_2020 <- t1.2027$`Basic payment scheme` + t1.2027$`Greening: practices beneficial for climate and environment`
#> Add new column showing 2027 payments
t1.2027 <- t1.2027 %>% 
   mutate(P1_2027 = case_when(P1_2020<=30000 ~ P1_2020* 0.15,
                             (P1_2020>30000 & P1_2020<=50000) ~ P1_2020* 0.15,
                             (P1_2020>50000 & P1_2020<=150000) ~ P1_2020* 0.10,
                             P1_2020> 150000 ~ P1_2020* 0.10))
# Group by postcode an summarise
t1.2027 <- t1.2027 %>% 
  group_by(PostcodePrefix_F202B) %>% 
  summarise(P1_2027 = sum(P1_2027))
#> Round
t1.2027$P1_2027 <- round(t1.2027$P1_2027, digits = 2)
# merge with main P1 results table (m.)
t1 <- merge(t1, t1.2027, by.x = "PostcodePrefix_F202B", by.y = "PostcodePrefix_F202B", all.x = TRUE)
# Change values according to agi land area proportion
t1$P1_2027 <- t1$P1_2027 * t1$PCENT_AGRI_LAND #> Round
t1$P1_2027 <- round(t1$P1_2027, digits = 2)
#> Check sum
sum.p1.2027 <- sum(t1$P1_2027)

#> Convert all NAs to zeros in data frame
t1[is.na(t1)] <- 0



#> Create plot

#> P1 payments by year
p1.in <- as.data.frame(colSums(t1[6:13]))
#> Change rownames to column "Value"
p1.plot <- tibble::rownames_to_column(p1.in, "VALUE")
#> Extract year from string
p1.plot$VALUE <- sub("^.*([0-9]{4}).*", "\\1", p1.plot$VALUE)
names(p1.plot)[1]<-paste("Year")
names(p1.plot)[2]<-paste("Value_bps")

#> Print table t1 on web page
print(knitr::kable(t1, col.names = c("PC", "2020_Unadjusted", "Agri_Total", "Agri_In", "Agri_PC", "2020", "2021", "2022", "2023", "2024", "2025", "2026", "2027"), caption = paste0("Table 1.", i , " ", "P1 payments by National Park - ",  active_park)))

#> Add line breaks
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")


# Simple bar chart
plot1 <- ggplot(p1.plot, aes(x=Year, y=Value_bps)) + 
  theme_bw() +
  geom_bar(stat = "identity", fill="#2c66b8", width = 0.6) +
  geom_text(aes(label= paste0(round(Value_bps / 1000000, digits = 1), " M")), vjust = 1.5, colour = "white") +
  labs(title = paste0("Pillar 1 payments:", " ", active_park), x = "Year", y = "Value (£)") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6))
plot1

#> Add line breaks
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")
cat("\n")

#> Render on web page
print(plot1)

#> Export the plot
ggsave(here("Out", "P1", "Nat_Park", "Plot", paste0("P1_payments_", active_park, ".png")))

#> Export final t1 table
write_csv(t1, here("Out", "P1", "Nat_Park", "CSV", paste0("P1_payments_", active_park, ".csv")))

}
Table 1.1 P1 payments by National Park - Exmoor
PC 2020_Unadjusted Agri_Total Agri_In Agri_PC 2020 2021 2022 2023 2024 2025 2026 2027
EX31 4437568.5 217.860 51.85 0.238 1056141.3 952290.6 793869.4 635448.17 477026.97 371412.84 238660.35 145733.71
EX32 2503991.1 115.188 9.89 0.086 215343.2 194261.1 161959.6 129658.11 97356.62 75822.30 48377.68 29572.42
EX34 1182067.8 74.954 9.28 0.124 146576.4 135698.4 113711.9 91725.47 69739.01 55081.36 35659.71 21494.27
EX35 653966.0 64.790 64.79 1.000 653966.0 574485.0 476390.1 378295.20 280200.30 214803.70 134980.84 83839.57
EX36 4654512.0 230.731 46.98 0.204 949520.4 858541.3 716113.2 573685.18 431257.11 336305.06 214910.34 131193.18
TA22 1853621.9 96.166 71.84 0.747 1384655.6 1234557.1 1026858.8 819160.46 611462.12 472996.57 306159.62 187696.20
TA23 740420.7 48.768 31.32 0.642 475350.1 438343.8 367041.3 295738.76 224436.24 176901.24 113860.48 68814.00
TA24 3715605.7 308.095 289.97 0.941 3496385.0 3091665.3 2567207.6 2042749.81 1518292.06 1168653.55 753364.72 464091.99
TA4 4151861.2 262.718 17.50 0.067 278174.7 253424.1 211697.9 169971.67 128245.46 100427.99 64168.87 39038.80

Table 1.2 P1 payments by National Park - Dartmoor
PC 2020_Unadjusted Agri_Total Agri_In Agri_PC 2020 2021 2022 2023 2024 2025 2026 2027
EX20 7462904.4 443.972 109.29 0.246 1835874.48 1663779.96 1388398.79 1113017.62 837636.45 654049.00 414743.79 253268.76
EX6 2976601.0 173.859 45.01 0.259 770939.65 693359.70 577718.75 462077.81 346436.86 269342.89 174719.42 106633.20
PL19 2872716.9 178.363 91.03 0.510 1465085.62 1304673.84 1084911.00 865148.15 645385.31 498876.75 318825.43 196039.86
PL20 2201678.1 241.356 207.84 0.861 1895644.83 1645430.69 1361083.97 1076737.25 792390.53 602826.04 385916.38 240349.31
PL21 1668995.0 118.524 46.83 0.395 659253.02 603604.19 504716.23 405828.28 306940.33 241015.02 155203.88 94083.26
PL6 358828.4 13.705 2.46 0.179 64230.28 53263.72 43629.18 33994.63 24360.09 17937.07 12030.83 7621.17
PL7 494378.8 47.727 18.25 0.382 188852.71 170610.06 142282.16 113954.25 85626.34 66741.07 42176.56 25809.60
TQ10 545741.0 57.013 36.93 0.648 353640.17 329140.84 276094.81 223048.79 170002.76 134638.74 84941.89 51311.95
TQ11 581560.0 49.024 37.52 0.765 444893.42 392806.30 326072.28 259338.27 192604.26 148114.92 94287.18 58265.93
TQ12 1454953.6 89.134 10.46 0.117 170229.57 156000.90 130466.46 104932.02 79397.59 62374.63 40328.73 24420.10
TQ13 4300642.1 302.455 238.90 0.790 3397507.25 3084325.87 2574699.78 2065073.69 1555447.61 1215696.88 775695.25 472785.30

4. Agri-Environment

4.1 ES - County-level

County-level analysis of Environmental Stewardship (ES) payments. Methods overview to follow.

#> Get ES data
#> This dataset has been extracted from the latest 914/02/2022) CS data for England
#> The pole of inaccessibility for each ES polygon has been pre-computed in QGIS, and the points then clipped to the extent of the SW LEP to reduce processing overheads
ES <- st_read(here("In", "Shape", "AES", "ES_SouthWest_PoI.shp"), stringsAsFactors = FALSE, quiet = TRUE)

i <- 0
#> Initiate loop
for(active_county in counties){
  
#> Loop counter (for table caption number)
i <- i+1

#> Get a boundary shapefile for the active county
target_county <- counties.sw |>
  filter(NAME == active_county)


# #> Test county polygon (Cornwall) - need to automate for all counties - then districts
# test.county <- cornwall.sf

#> Intersect with target county polygon
int <- st_intersection(ES, target_county)

#> Create non-geom version 
ES.ng <- int
st_geometry(ES.ng) <- NULL

#> Filter records to get ones that were live in 2020
#> New start date field in date format
ES.ng$STARTDAT.2 <- as.Date(ES.ng$STARTDAT, format = "%d/%m/%Y") # correct lowercase y
#> New start date field in date format
ES.ng$ENDDATE.2 <- as.Date(ES.ng$ENDDATE, format = "%d/%m/%Y") # correct lowercase y

#> Create column for length of agreement in years
ES.ng$Years <- as.numeric(ES.ng$ENDDATE.2 - ES.ng$STARTDAT.2) / 365
ES.ng$Years <- round(ES.ng$Years, 0)
ES.ng$TOTCOST <- as.numeric(ES.ng$TOTCOST)

#> Create column for average cost
ES.ng$AV.YEAR.COST <- ES.ng$TOTCOST / ES.ng$Years
#> Change NAs to Zero
ES.ng$AV.YEAR.COST[is.na(ES.ng$AV.YEAR.COST)] <- 0

#> Export as CSV
write_csv(ES.ng, here("Out", "AES", "ES", "County", paste0("ES_Agreements_", active_county, ".csv")))


#> 2018 and 2019 baseline
#> Calculate the "baseline" for agreements live in 2018 - earliest end date is "2020-04-30" - so all records will be included
live.2018 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(STARTDAT.2 <= "2017-12-31", ENDDATE.2 <= "2040-12-31")   
#> Sum Av cost for annual total for 2018
x.2018 <- sum(live.2018$AV.YEAR.COST)
x.2019 <- sum(live.2018$AV.YEAR.COST)  #2019 is the same - as same number of live agreements running throughout 2019-2018
#> Create new dataframe to hold this
totals <- c(x.2018, x.2019)
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0)))
#> Add figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)



#> Work out 2020 totals
#> Filter for agreements ending in 2020
end.2020 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(STARTDAT.2 <= "2017-12-31", ENDDATE.2 <= "2020-12-31")
#> Work out how many days in year that agreement running for
end.2020$YearStart <- as.Date("2020-01-01")
end.2020$Days <- as.numeric(end.2020$ENDDATE.2 - end.2020$YearStart)
#> Work out proportion of year agreement is active
end.2020$Year_Prop <- end.2020$Days / 365
#> New column for cost of 2020 based on proportion
end.2020$Cost_Actual <- end.2020$Year_Prop * end.2020$AV.YEAR.COST
#> Total amount for agreements ending in 2020
tot.2020 <- sum(end.2020$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2020 <- x.2018 - tot.2020
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)



#> Work out 2021 totals
#> Filter for agreements ending in 2021
end.2021 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2021-01-01", ENDDATE.2 <= "2021-12-31")
#> Work out how many days in year that agreement running for
end.2021$YearStart <- as.Date("2021-01-01")
end.2021$Days <- as.numeric(end.2021$ENDDATE.2 - end.2021$YearStart)
#> Work out proportion of year agreement is active
end.2021$Year_Prop <- end.2021$Days / 365
#> New column for cost of 2021 based on proportion
end.2021$Cost_Actual <- end.2021$Year_Prop * end.2021$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2021 <- sum(end.2021$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2021 <- x.2018 - (tot.2020 + tot.2021)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021)
#> Insert figure for 2021
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)



#> Work out 2022 totals
#> Filter for agreements ending in 2022
end.2022 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2022-01-01", ENDDATE.2 <= "2022-12-31")
#> Work out how mnay days in year that agreement running for
end.2022$YearStart <- as.Date("2022-01-01")
end.2022$Days <- as.numeric(end.2022$ENDDATE.2 - end.2022$YearStart)
#> Work out proipertion of year agreement is active
end.2022$Year_Prop <- end.2022$Days / 365
#> New column for cost of 2021 based on proportion
end.2022$Cost_Actual <- end.2022$Year_Prop * end.2022$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2022 <- sum(end.2022$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2022  <- x.2018 - (tot.2021 + tot.2020 + tot.2022)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)


#> Work out 2023 totals
#> Filter for agreements ending in 2023
end.2023 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2023-01-01", ENDDATE.2 <= "2023-12-31")
#> Work out how mnay days in year that agreement running for
end.2023$YearStart <- as.Date("2023-01-01")
end.2023$Days <- as.numeric(end.2023$ENDDATE.2 - end.2023$YearStart)
#> Work out proipertion of year agreement is active
end.2023$Year_Prop <- end.2023$Days / 365
#> New column for cost of 2021 based on proportion
end.2023$Cost_Actual <- end.2023$Year_Prop * end.2023$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2023 <- sum(end.2023$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2023  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)


#> Work out 2024 totals
#> Filter for agreements ending in 2024
end.2024 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2024-01-01", ENDDATE.2 <= "2024-12-31")
#> Work out how mnay days in year that agreement running for
end.2024$YearStart <- as.Date("2024-01-01")
end.2024$Days <- as.numeric(end.2024$ENDDATE.2 - end.2024$YearStart)
#> Work out proipertion of year agreement is active
end.2024$Year_Prop <- end.2024$Days / 365
#> New column for cost of 2021 based on proportion
end.2024$Cost_Actual <- end.2024$Year_Prop * end.2024$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2024 <- sum(end.2024$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2024  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023 + tot.2024)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0), Total_2024 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023, x.2024)
#> Insert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)


# Export ES payments summary table
write_csv(df, here("Out", "AES", "ES", "County", "CSV", paste0("ES_payments_Summary_", active_county, ".csv")))

#> Print table district.summary on web page
print(knitr::kable(df, caption = paste0("Table 4.1.", i , " ", "ES payments - ",  active_county)))


#> Plot data as bar chart
#> Make long
df.plot <- melt(df)
#> Extract year from string
names(df.plot)[1]<-paste("Year")
names(df.plot)[2]<-paste("Value_es")
df.plot$Year <- sub("^.*([0-9]{4}).*", "\\1", df.plot$Year)

# Simple bar chart
plot.es <- ggplot(df.plot, aes(x=Year, y=Value_es)) + 
  geom_bar(stat = "identity", fill="#b3a96fff", width = 0.6) +
  theme_bw() +
  labs(title = paste0("Environmental Stewardship payments: ", active_county), x = "Year", y = "Value (£)") +
  geom_text(aes(label= paste0(round(Value_es / 1000000, digits = 2), " M")), vjust = 1.5, colour = "white") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
plot.es
print(plot.es)
#> Export the plot
ggsave(here("Out", "AES", "ES", "County", "Plot", paste0("ES_payments_", active_county, ".png")), width = 7, height = 5)


#> Render table on web page



}
Table 4.1.1 ES payments - Cornwall
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
4712443 4712443 4016560 3562409 3082647 2147535 1846658

Table 4.1.2 ES payments - Devon
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
11013951 11013951 9952321 8453805 7227957 5902100 5056298

Table 4.1.3 ES payments - Dorset
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
3792714 3792714 3530843 2985831 2198648 1667019 1434688

Table 4.1.4 ES payments - Somerset
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
6850589 6850589 5914083 4953293 4382361 3639042 3261606

4.2 ES - District-level

District-level analysis of Environmental Stewardship (ES) payments. Methods overview to follow.

#> Get ES data
#> This dataset has been extracted from the latest 914/02/2022) CS data for England
#> The pole of inaccessibility for each ES polygon has been pre-computed in QGIS, and the points then clipped to the extent of the SW LEP to reduce processing overheads
ES <- st_read(here("In", "Shape", "AES", "ES_SouthWest_PoI.shp"), stringsAsFactors = FALSE, quiet = TRUE)

#> Get district boundary data for SW counties
districts.sw <- st_read(here("In", "Shape", "Districts.SW.UPDATED.shp"), quiet = TRUE)

#? Remove string "(B)" from district names
# # districts.sw$NAME <- str_remove_all(districts.sw$NAME, "(B)")
districts.sw$NAME<- gsub("\\s*\\([^\\)]+\\)", "", districts.sw$NAME)
#> Trim WS
districts.sw$NAME<- trimws(districts.sw$NAME)
# glimpse(districts.sw)
  
#> Create list for loop
districts.list <- as.list(districts.sw$NAME)

i <- 0
#> Initiate loop
for(active_district in districts.list){

#> Loop counter (for table caption number)
i <- i+1

#> Get a boundary shapefile for the active county
target_district <- districts.sw |>
  filter(NAME == active_district)


# #> Test county polygon (Cornwall) - need to automate for all counties - then districts
# test.county <- cornwall.sf

#> Intersect with target county polygon
int <- st_intersection(ES, target_district)

#> Create non-geom version 
ES.ng <- int
st_geometry(ES.ng) <- NULL

#> Filter records to get ones that were live in 2020
#> New start date field in date format
ES.ng$STARTDAT.2 <- as.Date(ES.ng$STARTDAT, format = "%d/%m/%Y") # correct lowercase y
#> New start date field in date format
ES.ng$ENDDATE.2 <- as.Date(ES.ng$ENDDATE, format = "%d/%m/%Y") # correct lowercase y

#> Create column for length of agreement in years
ES.ng$Years <- as.numeric(ES.ng$ENDDATE.2 - ES.ng$STARTDAT.2) / 365
ES.ng$Years <- round(ES.ng$Years, 0)
ES.ng$TOTCOST <- as.numeric(ES.ng$TOTCOST)

#> Create column for average cost
ES.ng$AV.YEAR.COST <- ES.ng$TOTCOST / ES.ng$Years
#> Change NAs to Zero
ES.ng$AV.YEAR.COST[is.na(ES.ng$AV.YEAR.COST)] <- 0

#> Export as CSV
write_csv(ES.ng, here("Out", "AES", "ES", "District", "CSV", paste0("ES_Agreements_", active_district, ".csv")))


#> 2018 and 2019 baseline
#> Calculate the "baseline" for agreements live in 2018 - earliest end date is "2020-04-30" - so all records will be included
live.2018 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(STARTDAT.2 <= "2017-12-31", ENDDATE.2 <= "2040-12-31")   
#> Sum Av cost for annual total for 2018
x.2018 <- sum(live.2018$AV.YEAR.COST)
x.2019 <- sum(live.2018$AV.YEAR.COST)  #2019 is the same - as same number of live agreements running throughout 2019-2018
#> Create new dataframe to hold this
totals <- c(x.2018, x.2019)
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0)))
#> Add figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)



#> Work out 2020 totals
#> Filter for agreements ending in 2020
end.2020 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(STARTDAT.2 <= "2017-12-31", ENDDATE.2 <= "2020-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2020)[1] != 0) {
#> Work out how many days in year that agreement running for
end.2020$YearStart <- as.Date("2020-01-01")
end.2020$Days <- as.numeric(end.2020$ENDDATE.2 - end.2020$YearStart)
#> Work out proportion of year agreement is active
end.2020$Year_Prop <- end.2020$Days / 365
#> New column for cost of 2020 based on proportion
end.2020$Cost_Actual <- end.2020$Year_Prop * end.2020$AV.YEAR.COST
#> Total amount for agreements ending in 2020
tot.2020 <- sum(end.2020$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2020 <- x.2018 - tot.2020
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

} else {
  
tot.2020 = 0
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2020 <- x.2018 - tot.2020
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

}





#> Work out 2021 totals
#> Filter for agreements ending in 2021
end.2021 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2021-01-01", ENDDATE.2 <= "2021-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2021)[1] != 0) {
#> Work out how many days in year that agreement running for
end.2021$YearStart <- as.Date("2021-01-01")
end.2021$Days <- as.numeric(end.2021$ENDDATE.2 - end.2021$YearStart)
#> Work out proportion of year agreement is active
end.2021$Year_Prop <- end.2021$Days / 365
#> New column for cost of 2021 based on proportion
end.2021$Cost_Actual <- end.2021$Year_Prop * end.2021$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2021 <- sum(end.2021$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2021 <- x.2018 - (tot.2020 + tot.2021)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021)
#> Insert figure for 2021
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

} else {
  
end.2021 = 0
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2021 <- x.2018 - (tot.2020 + tot.2021)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021)
#> Insert figure for 2021
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

}


#> Work out 2022 totals
#> Filter for agreements ending in 2022
end.2022 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2022-01-01", ENDDATE.2 <= "2022-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2022)[1] != 0) {
#> Work out how mnay days in year that agreement running for
end.2022$YearStart <- as.Date("2022-01-01")
end.2022$Days <- as.numeric(end.2022$ENDDATE.2 - end.2022$YearStart)
#> Work out proipertion of year agreement is active
end.2022$Year_Prop <- end.2022$Days / 365
#> New column for cost of 2021 based on proportion
end.2022$Cost_Actual <- end.2022$Year_Prop * end.2022$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2022 <- sum(end.2022$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2022  <- x.2018 - (tot.2021 + tot.2020 + tot.2022)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

} else {
end.2022 = 0
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2022  <- x.2018 - (tot.2021 + tot.2020 + tot.2022)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

}



#> Work out 2023 totals
#> Filter for agreements ending in 2023
end.2023 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2023-01-01", ENDDATE.2 <= "2023-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2023)[1] != 0) {
#> Work out how mnay days in year that agreement running for
end.2023$YearStart <- as.Date("2023-01-01")
end.2023$Days <- as.numeric(end.2023$ENDDATE.2 - end.2023$YearStart)
#> Work out proipertion of year agreement is active
end.2023$Year_Prop <- end.2023$Days / 365
#> New column for cost of 2021 based on proportion
end.2023$Cost_Actual <- end.2023$Year_Prop * end.2023$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2023 <- sum(end.2023$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2023  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)
} else {
  
end.2023 = 0
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2023  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)
}


#> Work out 2024 totals
#> Filter for agreements ending in 2024
end.2024 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2024-01-01", ENDDATE.2 <= "2024-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2024)[1] != 0) {
#> Work out how mnay days in year that agreement running for
end.2024$YearStart <- as.Date("2024-01-01")
end.2024$Days <- as.numeric(end.2024$ENDDATE.2 - end.2024$YearStart)
#> Work out proipertion of year agreement is active
end.2024$Year_Prop <- end.2024$Days / 365
#> New column for cost of 2021 based on proportion
end.2024$Cost_Actual <- end.2024$Year_Prop * end.2024$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2024 <- sum(end.2024$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2024  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023 + tot.2024)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0), Total_2024 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023, x.2024)
#> Insert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)
} else {
  
end.2024 = 0
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2024  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023 + tot.2024)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0), Total_2024 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023, x.2024)

}


# Export ES payments summary table
write_csv(df, here("Out", "AES", "ES", "District", "CSV", paste0("ES_payments_Summary_", active_district, ".csv")))

#> Print table district.summary on web page
print(knitr::kable(df, caption = paste0("Table 4.2.", i , " ", "ES payments - ",  active_district)))


#> Plot data as bar chart
#> Make long
df.plot <- melt(df)
#> Extract year from string
names(df.plot)[1]<-paste("Year")
names(df.plot)[2]<-paste("Value_es")
df.plot$Year <- sub("^.*([0-9]{4}).*", "\\1", df.plot$Year)

# Simple bar chart
plot.es <- ggplot(df.plot, aes(x=Year, y=Value_es)) + 
  geom_bar(stat = "identity", fill="#b3a96fff", width = 0.6) +
  theme_bw() +
  labs(title = paste0("Environmental Stewardship payments: ", active_district), x = "Year", y = "Value (£)") +
  geom_text(aes(label= paste0(round(Value_es / 1000000, digits = 2), " M")), vjust = 1.5, colour = "white") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
plot.es
print(plot.es)
#> Export the plot
ggsave(here("Out", "AES", "ES", "District", "Plot", paste0("ES_payments_", active_district, ".png")), width = 7, height = 5)




}
Table 4.2.1 ES payments - Bournemouth, Christchurch and Poole
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
106083.3 106083.3 106083.3 105745.8 95076.05 56370.75 42983.28

Table 4.2.2 ES payments - Dorset
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
3686631 3686631 3424760 2880085 2103572 1610649 1391705

Table 4.2.3 ES payments - Torbay
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023
130117.7 130117.7 130117.7 97288.55 -679225.1 -743023.9

Table 4.2.4 ES payments - City of Plymouth
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023
28814.88 28814.88 14525.85 -18303.34 -794817 -858615.8

Table 4.2.5 ES payments - North Somerset
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
359381.6 359381.6 322044.4 255110.7 217396.2 170338.2 124191.5

Table 4.2.6 ES payments - Cornwall
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
4493822 4493822 3797939 3350640 2884667 2094653 1795709

Table 4.2.7 ES payments - Isles of Scilly
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
218621.1 218621.1 218621.1 211768.9 197979.5 52882.1 50948.91

Table 4.2.8 ES payments - Mendip District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
713945 713945 629311.8 520118.5 385887.9 345577.6 240947.4

Table 4.2.9 ES payments - North Devon District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1540236 1540236 1282486 1117967 1061165 751356 666864.6

Table 4.2.10 ES payments - East Devon District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
841154.4 841154.4 705007.6 587631.8 552878.5 400127.7 353576.3

Table 4.2.11 ES payments - Teignbridge District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1294813 1294813 1099954 985458.2 861186.7 616574.6 480444.9

Table 4.2.12 ES payments - West Devon District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
4072922 4072922 3728240 3420375 2919589 2656549 2287807

Table 4.2.13 ES payments - Mid Devon District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
330223.1 330223.1 326400.8 204392.7 180357.5 158980.4 132963.4

Table 4.2.14 ES payments - Exeter District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023
13440.89 13440.89 5634.126 -116373.9 -140409.2 -161786.2

Table 4.2.15 ES payments - Somerset West and Taunton District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
3221566 3221566 2623982 2233932 2177044 1755655 1585541

Table 4.2.16 ES payments - South Somerset District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1009229 1009229 933445.1 742860.9 600144 479407.4 446463.5

Table 4.2.17 ES payments - Sedgemoor District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1546468 1546468 1405300 1201271 1001888 888063.8 864462.2

Table 4.2.18 ES payments - South Hams District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1895441 1895441 1826524 1250439 926890.9 771803.4 698716.2

Table 4.2.19 ES payments - Torridge District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
866788.5 866788.5 833430.4 770091.7 608440.8 493059.6 382276.5

4.3 ES - National Parks

National Park-level analysis of Environmental Stewardship (ES) payments. Methods overview to follow.

#> Get ES data
#> This dataset has been extracted from the latest 914/02/2022) CS data for England
#> The pole of inaccessibility for each ES polygon has been pre-computed in QGIS, and the points then clipped to the extent of the SW LEP to reduce processing overheads
ES <- st_read(here("In", "Shape", "AES", "ES_SouthWest_PoI.shp"), stringsAsFactors = FALSE, quiet = TRUE)

#> Create filter string (get names first!)
nat.parks.list <- c("Exmoor", "Dartmoor")

#> Sequential integer counter (for figure numbers)
i <-0 

#> Initiate loop
for(active_park in nat.parks.list){
  
#> Loop counter (for table caption number)
i <- i+1

#> Select target county (this will be start of loop)
target_park <- nat.parks|> 
  filter(NAME == active_park)

# #> Test county polygon (Cornwall) - need to automate for all counties - then districts
# test.county <- cornwall.sf

#> Intersect with target county polygon
int <- st_intersection(ES, target_park)

#> Create non-geom version 
ES.ng <- int
st_geometry(ES.ng) <- NULL

#> Filter records to get ones that were live in 2020
#> New start date field in date format
ES.ng$STARTDAT.2 <- as.Date(ES.ng$STARTDAT, format = "%d/%m/%Y") # correct lowercase y
#> New start date field in date format
ES.ng$ENDDATE.2 <- as.Date(ES.ng$ENDDATE, format = "%d/%m/%Y") # correct lowercase y

#> Create column for length of agreement in years
ES.ng$Years <- as.numeric(ES.ng$ENDDATE.2 - ES.ng$STARTDAT.2) / 365
ES.ng$Years <- round(ES.ng$Years, 0)
ES.ng$TOTCOST <- as.numeric(ES.ng$TOTCOST)

#> Create column for average cost
ES.ng$AV.YEAR.COST <- ES.ng$TOTCOST / ES.ng$Years
#> Change NAs to Zero
ES.ng$AV.YEAR.COST[is.na(ES.ng$AV.YEAR.COST)] <- 0

#> Export as CSV
write_csv(ES.ng, here("Out", "AES", "ES", "Nat_Park", "CSV", paste0("ES_Agreements_", active_park, ".csv")))


#> 2018 and 2019 baseline
#> Calculate the "baseline" for agreements live in 2018 - earliest end date is "2020-04-30" - so all records will be included
live.2018 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(STARTDAT.2 <= "2017-12-31", ENDDATE.2 <= "2040-12-31")   
#> Sum Av cost for annual total for 2018
x.2018 <- sum(live.2018$AV.YEAR.COST)
x.2019 <- sum(live.2018$AV.YEAR.COST)  #2019 is the same - as same number of live agreements running throughout 2019-2018
#> Create new dataframe to hold this
totals <- c(x.2018, x.2019)
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0)))
#> Add figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)



#> Work out 2020 totals
#> Filter for agreements ending in 2020
end.2020 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(STARTDAT.2 <= "2017-12-31", ENDDATE.2 <= "2020-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2020)[1] != 0) {
#> Work out how many days in year that agreement running for
end.2020$YearStart <- as.Date("2020-01-01")
end.2020$Days <- as.numeric(end.2020$ENDDATE.2 - end.2020$YearStart)
#> Work out proportion of year agreement is active
end.2020$Year_Prop <- end.2020$Days / 365
#> New column for cost of 2020 based on proportion
end.2020$Cost_Actual <- end.2020$Year_Prop * end.2020$AV.YEAR.COST
#> Total amount for agreements ending in 2020
tot.2020 <- sum(end.2020$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2020 <- x.2018 - tot.2020
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

} else {
  
tot.2020 = 0
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2020 <- x.2018 - tot.2020
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

}





#> Work out 2021 totals
#> Filter for agreements ending in 2021
end.2021 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2021-01-01", ENDDATE.2 <= "2021-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2021)[1] != 0) {
#> Work out how many days in year that agreement running for
end.2021$YearStart <- as.Date("2021-01-01")
end.2021$Days <- as.numeric(end.2021$ENDDATE.2 - end.2021$YearStart)
#> Work out proportion of year agreement is active
end.2021$Year_Prop <- end.2021$Days / 365
#> New column for cost of 2021 based on proportion
end.2021$Cost_Actual <- end.2021$Year_Prop * end.2021$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2021 <- sum(end.2021$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2021 <- x.2018 - (tot.2020 + tot.2021)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021)
#> Insert figure for 2021
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

} else {
  
end.2021 = 0
#> Subtract from 2018 baseline total to calculate final total for 2020
x.2021 <- x.2018 - (tot.2020 + tot.2021)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021)
#> Insert figure for 2021
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

}


#> Work out 2022 totals
#> Filter for agreements ending in 2022
end.2022 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2022-01-01", ENDDATE.2 <= "2022-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2022)[1] != 0) {
#> Work out how mnay days in year that agreement running for
end.2022$YearStart <- as.Date("2022-01-01")
end.2022$Days <- as.numeric(end.2022$ENDDATE.2 - end.2022$YearStart)
#> Work out proipertion of year agreement is active
end.2022$Year_Prop <- end.2022$Days / 365
#> New column for cost of 2021 based on proportion
end.2022$Cost_Actual <- end.2022$Year_Prop * end.2022$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2022 <- sum(end.2022$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2022  <- x.2018 - (tot.2021 + tot.2020 + tot.2022)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

} else {
end.2022 = 0
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2022  <- x.2018 - (tot.2021 + tot.2020 + tot.2022)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)

}



#> Work out 2023 totals
#> Filter for agreements ending in 2023
end.2023 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2023-01-01", ENDDATE.2 <= "2023-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2023)[1] != 0) {
#> Work out how mnay days in year that agreement running for
end.2023$YearStart <- as.Date("2023-01-01")
end.2023$Days <- as.numeric(end.2023$ENDDATE.2 - end.2023$YearStart)
#> Work out proipertion of year agreement is active
end.2023$Year_Prop <- end.2023$Days / 365
#> New column for cost of 2021 based on proportion
end.2023$Cost_Actual <- end.2023$Year_Prop * end.2023$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2023 <- sum(end.2023$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2023  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)
} else {
  
end.2023 = 0
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2023  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023)
#> AInsert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)
}


#> Work out 2024 totals
#> Filter for agreements ending in 2024
end.2024 <- ES.ng %>%
  select(AGREF, TOTCOST, Years, AV.YEAR.COST, STARTDAT.2, ENDDATE.2, Years) %>% 
  filter(ENDDATE.2  >= "2024-01-01", ENDDATE.2 <= "2024-12-31")
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(end.2024)[1] != 0) {
#> Work out how mnay days in year that agreement running for
end.2024$YearStart <- as.Date("2024-01-01")
end.2024$Days <- as.numeric(end.2024$ENDDATE.2 - end.2024$YearStart)
#> Work out proipertion of year agreement is active
end.2024$Year_Prop <- end.2024$Days / 365
#> New column for cost of 2021 based on proportion
end.2024$Cost_Actual <- end.2024$Year_Prop * end.2024$AV.YEAR.COST
#> Total amount for agreements ending in 2021
tot.2024 <- sum(end.2024$Cost_Actual)
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2024  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023 + tot.2024)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0), Total_2024 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023, x.2024)
#> Insert figure for 2018
df <- miscTools::insertRow(ES_Summary, 1, totals)
df <- as.data.frame(df)
} else {
  
end.2024 = 0
#> Subtract from 2018 baseline total to calculate final total for 2021
x.2024  <- x.2018 - (tot.2021 + tot.2020 + tot.2022 + tot.2023 + tot.2024)
#> Update df
ES_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0), Total_2024 = numeric(0)))
totals <- c(x.2018, x.2019, x.2020, x.2021, x.2022, x.2023, x.2024)

}


# Export ES payments summary table
write_csv(df, here("Out", "AES", "ES", "Nat_Park", "CSV", paste0("ES_payments_Summary_", active_park, ".csv")))

#> Print table district.summary on web page
print(knitr::kable(df, caption = paste0("Table 4.3.", i , " ", "ES payments - ",  active_park)))


#> Plot data as bar chart
#> Make long
df.plot <- melt(df)
#> Extract year from string
names(df.plot)[1]<-paste("Year")
names(df.plot)[2]<-paste("Value_es")
df.plot$Year <- sub("^.*([0-9]{4}).*", "\\1", df.plot$Year)

# Simple bar chart
plot.es <- ggplot(df.plot, aes(x=Year, y=Value_es)) + 
  geom_bar(stat = "identity", fill="#b3a96fff", width = 0.6) +
  theme_bw() +
  labs(title = paste0("Environmental Stewardship payments: ", active_park), x = "Year", y = "Value (£)") +
  geom_text(aes(label= paste0(round(Value_es / 1000000, digits = 2), " M")), vjust = 1.5, colour = "white") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
plot.es
print(plot.es)
#> Export the plot
ggsave(here("Out", "AES", "ES", "Nat_Park", "Plot", paste0("ES_payments_", active_park, ".png")), width = 7, height = 5)




}
Table 4.3.1 ES payments - Exmoor
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
3223260 3223260 2525247 2173958 2169447 1679060 1538239

Table 4.3.2 ES payments - Dartmoor
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
5005046 5005046 4644541 3913283 3326141 3005470 2571464

4.4 CS - County-level

County-level analysis of Countryside Stewardship (CS) payments. Methods overview to follow.

#> Get data
CS <- st_read(here("In", "Shape", "AES", "CS_SouthWest_PoI.shp"), stringsAsFactors = FALSE, quiet = TRUE)

#> Initiate counter (for figure numbers)
i <- 0
#> Initiate loop
for(active_county in counties){
#> Loop counter (for table caption number)
i <- i+1
#> Get a boundary shapefile for the active county
target_county <- counties.sw |>
filter(NAME == active_county)

# #> TESTING - one county at a time
# target_county <- counties.gb |>
#   filter(NAME == "Dorset")


#> Intersect with target county polygon
int <- st_intersection(CS, target_county)

#> Rename intersected data
CS.ng <- int

#> New start date field in date format
CS.ng$STARTDATE.2 <- lubridate::dmy_hms(CS.ng$STARTDATE) 
#> New end date field in date format
CS.ng$ENDDATE.2 <- lubridate::dmy_hms(CS.ng$ENDDATE) 
#> Format dates
CS.ng$STARTDATE.2 <- as.Date(CS.ng$STARTDATE.2, format = "%d/%m/%Y")
CS.ng$ENDDATE.2 <- as.Date(CS.ng$ENDDATE.2, format = "%d/%m/%Y")

#> Create a column to get number of years
CS.ng$Years <- as.numeric(CS.ng$ENDDATE.2 - CS.ng$STARTDATE.2)
CS.ng$Years <- CS.ng$Years /365
CS.ng$Years <- round(CS.ng$Years, 0)
#> List of agreement duration in years - unique durations
years.unique <- unique(CS.ng$Years)

# #> Export table for team use
# write_csv(CS.ng, here("Out", "AES", "CS", "County", "CSV", paste0("CS_Agreements_", active_county, ".csv")))


#> Export table for team use
write_csv(CS.ng, here("Out", "AES", "CS", "County", "CSV", paste0("CS_Agreements_", active_county, "Test.csv")))


#Create individual year columns to assign agreements to different years based on start date and duration - code for that to follow below
CS.ng$y.2016 <- 0
CS.ng$y.2017 <- 0 
CS.ng$y.2018 <- 0
CS.ng$y.2019 <- 0
CS.ng$y.2020 <- 0
CS.ng$y.2021 <- 0
CS.ng$y.2022 <- 0
CS.ng$y.2023 <- 0
CS.ng$y.2024 <- 0

#> Get unique combinations of agreements-years
unique.agree2 <- count(CS.ng, STARTDATE.2, Years) %>% 
  ungroup()
#> Export
write_csv(unique.agree2, (here("Out", "AES", "CS", "CS_Unique_Agree-Year_for_checking.csv")))


#> 1. Agreements starting in 2016 running for 10 years
cs.2016.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2016-01-01") & STARTDATE.2 <= as.Date("2016-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2016.10)[1] != 0) {
#> Update relevant cols
cs.2016.10$y.2016 = 1
cs.2016.10$y.2017 = 1
cs.2016.10$y.2018 = 1
cs.2016.10$y.2019 = 1
cs.2016.10$y.2020 = 1
cs.2016.10$y.2021 = 1
cs.2016.10$y.2022 = 1
cs.2016.10$y.2023 = 1
cs.2016.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2016[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2016
CS.ng$y.2017[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2017
CS.ng$y.2018[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2018
CS.ng$y.2019[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2019
CS.ng$y.2020[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2020
CS.ng$y.2021[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2021
CS.ng$y.2022[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2022
CS.ng$y.2023[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2023
CS.ng$y.2024[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2024
}


#> 2. Agreements starting in 2017 running for 10 years
cs.2017.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2017-01-01") & STARTDATE.2 <= as.Date("2017-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2017.10)[1] != 0) {
#> Update relevant cols
cs.2017.10$y.2017 = 1
cs.2017.10$y.2018 = 1
cs.2017.10$y.2019 = 1
cs.2017.10$y.2020 = 1
cs.2017.10$y.2021 = 1
cs.2017.10$y.2022 = 1
cs.2017.10$y.2023 = 1
cs.2017.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2017[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2017
CS.ng$y.2018[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2018
CS.ng$y.2019[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2019
CS.ng$y.2020[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2020
CS.ng$y.2021[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2021
CS.ng$y.2022[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2022
CS.ng$y.2023[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2023
CS.ng$y.2024[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2024
}



#> 3. Agreements starting in 2018 running for 5 years 
cs.2018.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2018-01-01") & STARTDATE.2 <= as.Date("2018-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2018.5)[1] != 0) {
#> Update relevant cols
cs.2018.5$y.2018 = 1
cs.2018.5$y.2019 = 1
cs.2018.5$y.2020 = 1
cs.2018.5$y.2021 = 1
cs.2018.5$y.2022 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2018[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2018
CS.ng$y.2019[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2019
CS.ng$y.2020[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2020
CS.ng$y.2021[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2021
CS.ng$y.2022[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2022
}



#> 4. Agreements starting in 2018 running for 10 years 
cs.2018.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2018-01-01") & STARTDATE.2 <= as.Date("2018-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2018.10)[1] != 0) {
#> Update relevant cols
cs.2018.10$y.2018 = 1
cs.2018.10$y.2019 = 1
cs.2018.10$y.2020 = 1
cs.2018.10$y.2021 = 1
cs.2018.10$y.2022 = 1
cs.2018.10$y.2023 = 1
cs.2018.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2018[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2018
CS.ng$y.2019[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2019
CS.ng$y.2020[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2020
CS.ng$y.2021[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2021
CS.ng$y.2022[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2022
CS.ng$y.2023[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2023
CS.ng$y.2024[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2024
}



#> 5. Agreements starting in 2019 running for 5 years 
cs.2019.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2019-01-01") & STARTDATE.2 <= as.Date("2019-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2019.5)[1] != 0) {
#> Update relevant cols
cs.2019.5$y.2019 = 1
cs.2019.5$y.2020 = 1
cs.2019.5$y.2021 = 1
cs.2019.5$y.2022 = 1
cs.2019.5$y.2023 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2019[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2019
CS.ng$y.2020[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2020
CS.ng$y.2021[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2021
CS.ng$y.2022[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2022
CS.ng$y.2023[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2023
}



#> 6. Agreements starting in 2019 running for 10 years 
cs.2019.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2019-01-01") & STARTDATE.2 <= as.Date("2019-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2019.10)[1] != 0) {
#> Update relevant cols
cs.2019.10$y.2019 = 1
cs.2019.10$y.2020 = 1
cs.2019.10$y.2021 = 1
cs.2019.10$y.2022 = 1
cs.2019.10$y.2023 = 1
cs.2019.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2019[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2019
CS.ng$y.2020[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2020
CS.ng$y.2021[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2021
CS.ng$y.2022[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2022
CS.ng$y.2023[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2023
CS.ng$y.2024[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2024
}



#> 7. Agreements starting in 20 running for 2years 
cs.2020.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.2)[1] != 0) {
#> Update relevant cols
cs.2020.2$y.2020 = 1
cs.2020.2$y.2021 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.2$CSREF,  CS.ng$CSREF)] <- cs.2020.2$y.2020
CS.ng$y.2021[match(cs.2020.2$CSREF,  CS.ng$CSREF)] <- cs.2020.2$y.2021
}



#> 8. Agreements starting in 20 running for 5 years 
cs.2020.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.5)[1] != 0) {
#> Update relevant cols
cs.2020.5$y.2020 = 1
cs.2020.5$y.2021 = 1
cs.2020.5$y.2022 = 1
cs.2020.5$y.2023 = 1
cs.2020.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2020
CS.ng$y.2021[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2021
CS.ng$y.2022[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2022
CS.ng$y.2023[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2023
CS.ng$y.2024[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2024
}



#> 9. Agreements starting in 20 running for 10 years 
cs.2020.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.10)[1] != 0) {
#> Update relevant cols
cs.2020.10$y.2020 = 1
cs.2020.10$y.2021 = 1
cs.2020.10$y.2022 = 1
cs.2020.10$y.2023 = 1
cs.2020.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2020
CS.ng$y.2021[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2021
CS.ng$y.2022[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2022
CS.ng$y.2023[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2023
CS.ng$y.2024[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2024
}



#> 10. Agreements starting in 21 running for 2 years 
cs.2021.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.2)[1] != 0) {
#> Update relevant cols
cs.2021.2$y.2021 = 1
cs.2021.2$y.2022 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.2$CSREF,  CS.ng$CSREF)] <- cs.2021.2$y.2021
CS.ng$y.2022[match(cs.2021.2$CSREF,  CS.ng$CSREF)] <- cs.2021.2$y.2022
}


#> 11. Agreements starting in 21 running for 5 years 
cs.2021.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.5 )[1] != 0) {
#> Update relevant cols
cs.2021.5$y.2021 = 1
cs.2021.5$y.2022 = 1
cs.2021.5$y.2023 = 1
cs.2021.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2021
CS.ng$y.2022[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2022
CS.ng$y.2023[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2023
CS.ng$y.2024[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2024
}



#> 11. Agreements starting in 21 running for 10 years 
cs.2021.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.10)[1] != 0) {
#> Update relevant cols
cs.2021.10$y.2021 = 1
cs.2021.10$y.2022 = 1
cs.2021.10$y.2023 = 1
cs.2021.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2021
CS.ng$y.2022[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2022
CS.ng$y.2023[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2023
CS.ng$y.2024[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2024
}



#> 12. Agreements starting in 21 running for 20 years 
cs.2021.20 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 20)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.20 )[1] != 0) {
#> Update relevant cols
cs.2021.20$y.2021 = 1
cs.2021.20$y.2022 = 1
cs.2021.20$y.2023 = 1
cs.2021.20$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2021
CS.ng$y.2022[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2022
CS.ng$y.2023[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2023
CS.ng$y.2024[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2024
}



#> 13. Agreements starting in 22 running for 2 years 
cs.2022.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.2)[1] != 0) {
#> Update relevant cols
cs.2022.2$y.2022 = 1
cs.2022.2$y.2023 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.2$CSREF,  CS.ng$CSREF)] <- cs.2022.2$y.2022
CS.ng$y.2023[match(cs.2022.2$CSREF,  CS.ng$CSREF)] <- cs.2022.2$y.2023
}



#> 14. Agreements starting in 22 running for 5 years 
cs.2022.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.5)[1] != 0) {
#> Update relevant cols
cs.2022.5$y.2022 = 1
cs.2022.5$y.2023 = 1
cs.2022.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2022
CS.ng$y.2023[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2023
CS.ng$y.2024[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2024
}



#> 15. Agreements starting in 22 running for 10 years 
cs.2022.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.10)[1] != 0) {
#> Update relevant cols
cs.2022.10$y.2022 = 1
cs.2022.10$y.2023 = 1
cs.2022.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2022
CS.ng$y.2023[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2023
CS.ng$y.2024[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2024
}



# Remove records where average annual cost id "unavailable"
CS.final <- CS.ng[!(CS.ng$AVGANNCOST=="*Unavailable"),]



# Export final raw data table
write_csv(CS.final, here("Out", "AES", "CS", "County", "CSV", paste0("CS_", active_county, "_Raw_Processed.csv")))

# # Export final raw data table - TESTING
# write_csv(CS.final, here("Out", "AES", "CS", "County", "CSV", paste0("CS__Raw_Processed_Test.csv")))



# Change annual cost column to numerical
CS.final$AVGANNCOST <- as.numeric(CS.final$AVGANNCOST)
CS.final$AVGANNCOST <- round(CS.final$AVGANNCOST, 2)


#> Calculate annual costs
#> Total for 2018
CS.2018.Total <- CS.final %>% 
  filter(y.2018 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2018 <- sum(CS.2018.Total$AVGANNCOST)

#> Total for 2019
CS.2019.Total <- CS.final %>% 
  filter(y.2019 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2019 <- sum(CS.2019.Total$AVGANNCOST)


#> Total for 2020
CS.2020.Total <- CS.final %>% 
  filter(y.2020 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2020 <- sum(CS.2020.Total$AVGANNCOST)


#> Total for 2021
CS.2021.Total <- CS.final %>% 
  filter(y.2021 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2021 <- sum(CS.2021.Total$AVGANNCOST)

#> Total for 2022
CS.2022.Total <- CS.final %>% 
  filter(y.2022 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2022 <- sum(CS.2022.Total$AVGANNCOST)

#> Total for 2023
CS.2023.Total <- CS.final %>% 
  filter(y.2023 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2023 <- sum(CS.2023.Total$AVGANNCOST)

#> Total for 2024
CS.2024.Total <- CS.final %>% 
  filter(y.2024 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2024 <- sum(CS.2024.Total$AVGANNCOST)


#> Create data frame of annual costs
CS_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0), Total_2024 = numeric(0)))
totals.cs <- c(CS.sum.2018, CS.sum.2019, CS.sum.2020, CS.sum.2021, CS.sum.2022, CS.sum.2023, CS.sum.2024)
#> Insert figure for 2018
df2 <- miscTools::insertRow(CS_Summary, 1, totals.cs)
df2 <- as.data.frame(df2)

#> Print table district.summary on web page
print(knitr::kable(df2, caption = paste0("Table 4.4.", i , " ", "CS payments - ",  active_county)))

#> Make long
df2.plot <- melt(df2)
#> Extract year from string
names(df2.plot)[1]<-paste("Year")
names(df2.plot)[2]<-paste("Value_cs")
df2.plot$Year <- sub("^.*([0-9]{4}).*", "\\1", df2.plot$Year)

# Simple bar chart
plot.cs <- ggplot(df2.plot, aes(x=Year, y=Value_cs)) + 
  geom_bar(stat = "identity", fill="#e38d0b", width = 0.6) +
  theme_bw() +
  labs(title = paste0("Countryside Stewardship payments: ", active_county), x = "Year", y = "Value (£)") +
  geom_text(aes(label= paste0(round(Value_cs / 1000000, digits = 2), " M")), vjust = 1.5, colour = "white", size = 3.5) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
plot.cs
print(plot.cs)


#> Export the plot
ggsave(here("Out", "AES", "CS", "County", "Plot", paste0("CS_payments_", active_county, ".png")), width = 7, height = 5)

# #> Export the plot - TESTING
# ggsave(here("Out", "AES", "CS", "County", "Plot", paste0("CS_payments_County_Test.png")), width = 7, height = 5)


#Export table
write_csv(df2, here("Out", "AES", "CS", "County", "CSV", paste0("CS_payments_", active_county, ".csv")))


# #Export table - TESTING
# write_csv(df2, here("Out", "AES", "CS", "County", "CSV", "CS_payments_County_Test.csv"))

}
Table 4.4.1 CS payments - Cornwall
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1642619 3815312 6052031 10541488 10474307 8346188 6205684

Table 4.4.2 CS payments - Devon
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
3903397 8948297 15208232 24477806 23782115 18413622 13582947

Table 4.4.3 CS payments - Dorset
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1791615 2971533 4256189 6984160 6976659 5333381 4335423

Table 4.4.4 CS payments - Somerset
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1368900 3493153 6045385 10232878 10089629 7627047 5671227

4.5 CS - District-level

District-level analysis of Countryside Stewardship (CS) payments. Methods overview to follow.

#> Get data
CS <- st_read(here("In", "Shape", "AES", "CS_SouthWest_PoI.shp"), stringsAsFactors = FALSE, quiet = TRUE)

#> Get district boundary data for SW counties
districts.sw <- st_read(here("In", "Shape", "Districts.SW.UPDATED.shp"), quiet = TRUE)

#? Remove string "(B)" from district names
# # districts.sw$NAME <- str_remove_all(districts.sw$NAME, "(B)")
districts.sw$NAME<- gsub("\\s*\\([^\\)]+\\)", "", districts.sw$NAME)
#> Trim WS
districts.sw$NAME<- trimws(districts.sw$NAME)
# glimpse(districts.sw)
  
#> Create list for loop
districts.list <- as.list(districts.sw$NAME)



i <- 0
#> Initiate loop
for(active_district in districts.list){

#> Loop counter (for table caption number)
i <- i+1

#> Get a boundary shapefile for the active county
target_district <- districts.sw |>
  filter(NAME == active_district)


#> Intersect with target county polygon
int <- st_intersection(CS, target_district)

#> Rename intersected data
CS.ng <- int

#> New start date field in date format
CS.ng$STARTDATE.2 <- lubridate::dmy_hms(CS.ng$STARTDATE) 
#> New end date field in date format
CS.ng$ENDDATE.2 <- lubridate::dmy_hms(CS.ng$ENDDATE) 
#> Format dates
CS.ng$STARTDATE.2 <- as.Date(CS.ng$STARTDATE.2, format = "%d/%m/%Y")
CS.ng$ENDDATE.2 <- as.Date(CS.ng$ENDDATE.2, format = "%d/%m/%Y")

#> Create a column to get number of years
CS.ng$Years <- as.numeric(CS.ng$ENDDATE.2 - CS.ng$STARTDATE.2)
CS.ng$Years <- CS.ng$Years /365
CS.ng$Years <- round(CS.ng$Years, 0)
#> List of agreement duration in years - unique duration
years.unique <- unique(CS.ng$Years)

# #> Export table for team use
# write_csv(CS.ng, here("Out", "AES", "CS", "County", "CSV", paste0("CS_Agreements_", active_county, ".csv")))


#> Export table for team use
write_csv(CS.ng, here("Out", "AES", "CS", "District", "CSV", paste0("CS_Agreements", active_district, ".csv")))


#Create individual year columns to assign agreements to different years based on start date and duration - code for that to follow below
CS.ng$y.2016 <- 0
CS.ng$y.2017 <- 0 
CS.ng$y.2018 <- 0
CS.ng$y.2019 <- 0
CS.ng$y.2020 <- 0
CS.ng$y.2021 <- 0
CS.ng$y.2022 <- 0
CS.ng$y.2023 <- 0
CS.ng$y.2024 <- 0

#> Get unique combinations of agreements-years
unique.agree2 <- count(CS.ng, STARTDATE.2, Years) %>% 
  ungroup()
#> Export
write_csv(unique.agree2, (here("Out", "AES", "CS", "CS_Unique_Agree-Year_for_checking.csv")))


#> 1. Agreements starting in 2016 running for 10 years
cs.2016.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2016-01-01") & STARTDATE.2 <= as.Date("2016-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2016.10)[1] != 0) {
#> Update relevant cols
cs.2016.10$y.2016 = 1
cs.2016.10$y.2017 = 1
cs.2016.10$y.2018 = 1
cs.2016.10$y.2019 = 1
cs.2016.10$y.2020 = 1
cs.2016.10$y.2021 = 1
cs.2016.10$y.2022 = 1
cs.2016.10$y.2023 = 1
cs.2016.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2016[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2016
CS.ng$y.2017[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2017
CS.ng$y.2018[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2018
CS.ng$y.2019[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2019
CS.ng$y.2020[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2020
CS.ng$y.2021[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2021
CS.ng$y.2022[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2022
CS.ng$y.2023[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2023
CS.ng$y.2024[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2024
}


#> 2. Agreements starting in 2017 running for 10 years
cs.2017.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2017-01-01") & STARTDATE.2 <= as.Date("2017-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2017.10)[1] != 0) {
#> Update relevant cols
cs.2017.10$y.2017 = 1
cs.2017.10$y.2018 = 1
cs.2017.10$y.2019 = 1
cs.2017.10$y.2020 = 1
cs.2017.10$y.2021 = 1
cs.2017.10$y.2022 = 1
cs.2017.10$y.2023 = 1
cs.2017.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2017[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2017
CS.ng$y.2018[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2018
CS.ng$y.2019[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2019
CS.ng$y.2020[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2020
CS.ng$y.2021[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2021
CS.ng$y.2022[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2022
CS.ng$y.2023[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2023
CS.ng$y.2024[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2024
}



#> 3. Agreements starting in 2018 running for 5 years 
cs.2018.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2018-01-01") & STARTDATE.2 <= as.Date("2018-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2018.5)[1] != 0) {
#> Update relevant cols
cs.2018.5$y.2018 = 1
cs.2018.5$y.2019 = 1
cs.2018.5$y.2020 = 1
cs.2018.5$y.2021 = 1
cs.2018.5$y.2022 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2018[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2018
CS.ng$y.2019[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2019
CS.ng$y.2020[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2020
CS.ng$y.2021[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2021
CS.ng$y.2022[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2022
}



#> 4. Agreements starting in 2018 running for 10 years 
cs.2018.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2018-01-01") & STARTDATE.2 <= as.Date("2018-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2018.10)[1] != 0) {
#> Update relevant cols
cs.2018.10$y.2018 = 1
cs.2018.10$y.2019 = 1
cs.2018.10$y.2020 = 1
cs.2018.10$y.2021 = 1
cs.2018.10$y.2022 = 1
cs.2018.10$y.2023 = 1
cs.2018.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2018[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2018
CS.ng$y.2019[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2019
CS.ng$y.2020[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2020
CS.ng$y.2021[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2021
CS.ng$y.2022[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2022
CS.ng$y.2023[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2023
CS.ng$y.2024[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2024
}



#> 5. Agreements starting in 2019 running for 5 years 
cs.2019.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2019-01-01") & STARTDATE.2 <= as.Date("2019-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2019.5)[1] != 0) {
#> Update relevant cols
cs.2019.5$y.2019 = 1
cs.2019.5$y.2020 = 1
cs.2019.5$y.2021 = 1
cs.2019.5$y.2022 = 1
cs.2019.5$y.2023 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2019[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2019
CS.ng$y.2020[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2020
CS.ng$y.2021[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2021
CS.ng$y.2022[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2022
CS.ng$y.2023[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2023
}



#> 6. Agreements starting in 2019 running for 10 years 
cs.2019.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2019-01-01") & STARTDATE.2 <= as.Date("2019-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2019.10)[1] != 0) {
#> Update relevant cols
cs.2019.10$y.2019 = 1
cs.2019.10$y.2020 = 1
cs.2019.10$y.2021 = 1
cs.2019.10$y.2022 = 1
cs.2019.10$y.2023 = 1
cs.2019.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2019[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2019
CS.ng$y.2020[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2020
CS.ng$y.2021[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2021
CS.ng$y.2022[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2022
CS.ng$y.2023[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2023
CS.ng$y.2024[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2024
}



#> 7. Agreements starting in 20 running for 2years 
cs.2020.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.2)[1] != 0) {
#> Update relevant cols
cs.2020.2$y.2020 = 1
cs.2020.2$y.2021 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.2$CSREF,  CS.ng$CSREF)] <- cs.2020.2$y.2020
CS.ng$y.2021[match(cs.2020.2$CSREF,  CS.ng$CSREF)] <- cs.2020.2$y.2021
}



#> 8. Agreements starting in 20 running for 5 years 
cs.2020.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.5)[1] != 0) {
#> Update relevant cols
cs.2020.5$y.2020 = 1
cs.2020.5$y.2021 = 1
cs.2020.5$y.2022 = 1
cs.2020.5$y.2023 = 1
cs.2020.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2020
CS.ng$y.2021[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2021
CS.ng$y.2022[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2022
CS.ng$y.2023[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2023
CS.ng$y.2024[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2024
}



#> 9. Agreements starting in 20 running for 10 years 
cs.2020.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.10)[1] != 0) {
#> Update relevant cols
cs.2020.10$y.2020 = 1
cs.2020.10$y.2021 = 1
cs.2020.10$y.2022 = 1
cs.2020.10$y.2023 = 1
cs.2020.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2020
CS.ng$y.2021[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2021
CS.ng$y.2022[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2022
CS.ng$y.2023[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2023
CS.ng$y.2024[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2024
}



#> 10. Agreements starting in 21 running for 2 years 
cs.2021.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.2)[1] != 0) {
#> Update relevant cols
cs.2021.2$y.2021 = 1
cs.2021.2$y.2022 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.2$CSREF,  CS.ng$CSREF)] <- cs.2021.2$y.2021
CS.ng$y.2022[match(cs.2021.2$CSREF,  CS.ng$CSREF)] <- cs.2021.2$y.2022
}


#> 11. Agreements starting in 21 running for 5 years 
cs.2021.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.5 )[1] != 0) {
#> Update relevant cols
cs.2021.5$y.2021 = 1
cs.2021.5$y.2022 = 1
cs.2021.5$y.2023 = 1
cs.2021.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2021
CS.ng$y.2022[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2022
CS.ng$y.2023[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2023
CS.ng$y.2024[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2024
}



#> 11. Agreements starting in 21 running for 10 years 
cs.2021.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.10)[1] != 0) {
#> Update relevant cols
cs.2021.10$y.2021 = 1
cs.2021.10$y.2022 = 1
cs.2021.10$y.2023 = 1
cs.2021.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2021
CS.ng$y.2022[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2022
CS.ng$y.2023[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2023
CS.ng$y.2024[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2024
}



#> 12. Agreements starting in 21 running for 20 years 
cs.2021.20 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 20)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.20 )[1] != 0) {
#> Update relevant cols
cs.2021.20$y.2021 = 1
cs.2021.20$y.2022 = 1
cs.2021.20$y.2023 = 1
cs.2021.20$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2021
CS.ng$y.2022[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2022
CS.ng$y.2023[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2023
CS.ng$y.2024[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2024
}



#> 13. Agreements starting in 22 running for 2 years 
cs.2022.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.2)[1] != 0) {
#> Update relevant cols
cs.2022.2$y.2022 = 1
cs.2022.2$y.2023 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.2$CSREF,  CS.ng$CSREF)] <- cs.2022.2$y.2022
CS.ng$y.2023[match(cs.2022.2$CSREF,  CS.ng$CSREF)] <- cs.2022.2$y.2023
}



#> 14. Agreements starting in 22 running for 5 years 
cs.2022.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.5)[1] != 0) {
#> Update relevant cols
cs.2022.5$y.2022 = 1
cs.2022.5$y.2023 = 1
cs.2022.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2022
CS.ng$y.2023[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2023
CS.ng$y.2024[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2024
}



#> 15. Agreements starting in 22 running for 10 years 
cs.2022.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.10)[1] != 0) {
#> Update relevant cols
cs.2022.10$y.2022 = 1
cs.2022.10$y.2023 = 1
cs.2022.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2022
CS.ng$y.2023[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2023
CS.ng$y.2024[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2024
}



# Remove records where average annual cost id "unavailable"
CS.final <- CS.ng[!(CS.ng$AVGANNCOST=="*Unavailable"),]



# Export final raw data table
write_csv(CS.final, here("Out", "AES", "CS", "District", "CSV", paste0("CS_", active_district, "_Raw_Processed.csv")))

# # Export final raw data table - TESTING
# write_csv(CS.final, here("Out", "AES", "CS", "County", "CSV", paste0("CS__Raw_Processed_Test.csv")))



# Change annual cost column to numerical
CS.final$AVGANNCOST <- as.numeric(CS.final$AVGANNCOST)
CS.final$AVGANNCOST <- round(CS.final$AVGANNCOST, 2)


#> Calculate annual costs
#> Total for 2018
CS.2018.Total <- CS.final %>% 
  filter(y.2018 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2018 <- sum(CS.2018.Total$AVGANNCOST)

#> Total for 2019
CS.2019.Total <- CS.final %>% 
  filter(y.2019 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2019 <- sum(CS.2019.Total$AVGANNCOST)


#> Total for 2020
CS.2020.Total <- CS.final %>% 
  filter(y.2020 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2020 <- sum(CS.2020.Total$AVGANNCOST)


#> Total for 2021
CS.2021.Total <- CS.final %>% 
  filter(y.2021 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2021 <- sum(CS.2021.Total$AVGANNCOST)

#> Total for 2022
CS.2022.Total <- CS.final %>% 
  filter(y.2022 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2022 <- sum(CS.2022.Total$AVGANNCOST)

#> Total for 2023
CS.2023.Total <- CS.final %>% 
  filter(y.2023 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2023 <- sum(CS.2023.Total$AVGANNCOST)

#> Total for 2024
CS.2024.Total <- CS.final %>% 
  filter(y.2024 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2024 <- sum(CS.2024.Total$AVGANNCOST)


#> Create data frame of annual costs
CS_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0), Total_2024 = numeric(0)))
totals.cs <- c(CS.sum.2018, CS.sum.2019, CS.sum.2020, CS.sum.2021, CS.sum.2022, CS.sum.2023, CS.sum.2024)
#> Insert figure for 2018
df2 <- miscTools::insertRow(CS_Summary, 1, totals.cs)
df2 <- as.data.frame(df2)

#> Print table district.summary on web page
print(knitr::kable(df2, caption = paste0("Table 4.5.", i , " ", "CS payments - ",  active_district)))

#> Make long
df2.plot <- melt(df2)
#> Extract year from string
names(df2.plot)[1]<-paste("Year")
names(df2.plot)[2]<-paste("Value_cs")
df2.plot$Year <- sub("^.*([0-9]{4}).*", "\\1", df2.plot$Year)

# Simple bar chart
plot.cs <- ggplot(df2.plot, aes(x=Year, y=Value_cs)) + 
  geom_bar(stat = "identity", fill="#e38d0b", width = 0.6) +
  theme_bw() +
  labs(title = paste0("Countryside Stewardship payments: ", active_district), x = "Year", y = "Value (£)") +
  geom_text(aes(label= paste0(round(Value_cs / 1000000, digits = 2), " M")), vjust = 1.5, colour = "white", size = 3.5) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
plot.cs
print(plot.cs)


#> Export the plot
ggsave(here("Out", "AES", "CS", "District", "Plot", paste0("CS_payments_", active_district, ".png")), width = 7, height = 5)

# #> Export the plot - TESTING
# ggsave(here("Out", "AES", "CS", "County", "Plot", paste0("CS_payments_County_Test.png")), width = 7, height = 5)


#Export table
write_csv(df2, here("Out", "AES", "CS", "District", "CSV", paste0("CS_payments_", active_district, ".csv")))


# #Export table - TESTING
# write_csv(df2, here("Out", "AES", "CS", "County", "CSV", "CS_payments_County_Test.csv"))

}
Table 4.5.1 CS payments - Bournemouth, Christchurch and Poole
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
103263.6 158993.2 194506.5 194506.5 194506.5 191512.7 135783

Table 4.5.2 CS payments - Dorset
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1688351 2812540 4061682 6789654 6782153 5141868 4199640

Table 4.5.3 CS payments - Torbay
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
0 0 18736.41 18736.41 18736.41 18736.41 18736.41

Table 4.5.4 CS payments - City of Plymouth
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
0 0 6242.79 6242.79 6242.79 6242.79 6242.79

Table 4.5.5 CS payments - North Somerset
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
39415.63 128626.8 223173.1 366769.6 368148.4 225785.6 175487.1

Table 4.5.6 CS payments - Cornwall
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1642619 3815312 6052031 10541488 10474307 8346188 6205684

Table 4.5.7 CS payments - Isles of Scilly
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
0 0 0 0 0 0 0

Table 4.5.8 CS payments - Mendip District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
180558 370796.8 681652.6 1417915 1411463 1066339 867417.9

Table 4.5.9 CS payments - North Devon District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
1054832 2186363 3714008 5570941 5378101 4017822 3002549

Table 4.5.10 CS payments - East Devon District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
624322.9 1593875 2536338 3879551 3803659 2955823 2016972

Table 4.5.11 CS payments - Teignbridge District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
183268.5 483544.2 890069.2 1613032 1554846 1216213 924866.4

Table 4.5.12 CS payments - West Devon District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
436003 1070140 1831753 3074129 2899122 2296230 1666072

Table 4.5.13 CS payments - Mid Devon District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
433334.3 1142100 1928984 3405391 3351289 2672281 1966826

Table 4.5.14 CS payments - Exeter District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
4105.02 22884.87 35852.96 35852.96 35852.96 31747.94 12968.09

Table 4.5.15 CS payments - Somerset West and Taunton District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
343245.8 1219330 2360442 3972653 3851860 3251261 2496645

Table 4.5.16 CS payments - South Somerset District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
641285.7 1314698 1945032 3051075 3001708 2057381 1409320

Table 4.5.17 CS payments - Sedgemoor District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
164395 459700.8 835084.9 1424465 1456451 1026281 722357.2

Table 4.5.18 CS payments - South Hams District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
452962.3 998096.5 1993245 2918487 2828357 2275218 1739957

Table 4.5.19 CS payments - Torridge District
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
714568.7 1451294 2253003 3955443 3905908 2923308 2227756

4.6 CS - National Parks

National Park-level analysis of Countryside Stewardship (CS) payments. Methods overview to follow.

#> Get data
CS <- st_read(here("In", "Shape", "AES", "CS_SouthWest_PoI.shp"), stringsAsFactors = FALSE, quiet = TRUE)

#> Create filter string (get names first!)
nat.parks.list <- c("Exmoor", "Dartmoor")

#> Sequential integer counter (for figure numbers)
i <-0 

#> Initiate loop
for(active_park in nat.parks.list){
  
#> Loop counter (for table caption number)
i <- i+1

#> Select target county (this will be start of loop)
target_park <- nat.parks|> 
  filter(NAME == active_park)


#> Intersect with target county polygon
int <- st_intersection(CS, target_park)

#> Rename intersected data
CS.ng <- int

#> New start date field in date format
CS.ng$STARTDATE.2 <- lubridate::dmy_hms(CS.ng$STARTDATE) 
#> New end date field in date format
CS.ng$ENDDATE.2 <- lubridate::dmy_hms(CS.ng$ENDDATE) 
#> Format dates
CS.ng$STARTDATE.2 <- as.Date(CS.ng$STARTDATE.2, format = "%d/%m/%Y")
CS.ng$ENDDATE.2 <- as.Date(CS.ng$ENDDATE.2, format = "%d/%m/%Y")

#> Create a column to get number of years
CS.ng$Years <- as.numeric(CS.ng$ENDDATE.2 - CS.ng$STARTDATE.2)
CS.ng$Years <- CS.ng$Years /365
CS.ng$Years <- round(CS.ng$Years, 0)
#> List of agreement duration in years - unique duration
years.unique <- unique(CS.ng$Years)

# #> Export table for team use
# write_csv(CS.ng, here("Out", "AES", "CS", "County", "CSV", paste0("CS_Agreements_", active_county, ".csv")))


#> Export table for team use
write_csv(CS.ng, here("Out", "AES", "CS", "Nat_Park", "CSV", paste0("CS_Agreements", active_park, ".csv")))


#Create individual year columns to assign agreements to different years based on start date and duration - code for that to follow below
CS.ng$y.2016 <- 0
CS.ng$y.2017 <- 0 
CS.ng$y.2018 <- 0
CS.ng$y.2019 <- 0
CS.ng$y.2020 <- 0
CS.ng$y.2021 <- 0
CS.ng$y.2022 <- 0
CS.ng$y.2023 <- 0
CS.ng$y.2024 <- 0

#> Get unique combinations of agreements-years
unique.agree2 <- count(CS.ng, STARTDATE.2, Years) %>% 
  ungroup()
#> Export
write_csv(unique.agree2, (here("Out", "AES", "CS", "CS_Unique_Agree-Year_for_checking.csv")))


#> 1. Agreements starting in 2016 running for 10 years
cs.2016.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2016-01-01") & STARTDATE.2 <= as.Date("2016-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2016.10)[1] != 0) {
#> Update relevant cols
cs.2016.10$y.2016 = 1
cs.2016.10$y.2017 = 1
cs.2016.10$y.2018 = 1
cs.2016.10$y.2019 = 1
cs.2016.10$y.2020 = 1
cs.2016.10$y.2021 = 1
cs.2016.10$y.2022 = 1
cs.2016.10$y.2023 = 1
cs.2016.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2016[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2016
CS.ng$y.2017[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2017
CS.ng$y.2018[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2018
CS.ng$y.2019[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2019
CS.ng$y.2020[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2020
CS.ng$y.2021[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2021
CS.ng$y.2022[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2022
CS.ng$y.2023[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2023
CS.ng$y.2024[match(cs.2016.10$CSREF,  CS.ng$CSREF)] <- cs.2016.10$y.2024
}


#> 2. Agreements starting in 2017 running for 10 years
cs.2017.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2017-01-01") & STARTDATE.2 <= as.Date("2017-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2017.10)[1] != 0) {
#> Update relevant cols
cs.2017.10$y.2017 = 1
cs.2017.10$y.2018 = 1
cs.2017.10$y.2019 = 1
cs.2017.10$y.2020 = 1
cs.2017.10$y.2021 = 1
cs.2017.10$y.2022 = 1
cs.2017.10$y.2023 = 1
cs.2017.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2017[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2017
CS.ng$y.2018[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2018
CS.ng$y.2019[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2019
CS.ng$y.2020[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2020
CS.ng$y.2021[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2021
CS.ng$y.2022[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2022
CS.ng$y.2023[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2023
CS.ng$y.2024[match(cs.2017.10$CSREF,  CS.ng$CSREF)] <- cs.2017.10$y.2024
}



#> 3. Agreements starting in 2018 running for 5 years 
cs.2018.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2018-01-01") & STARTDATE.2 <= as.Date("2018-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2018.5)[1] != 0) {
#> Update relevant cols
cs.2018.5$y.2018 = 1
cs.2018.5$y.2019 = 1
cs.2018.5$y.2020 = 1
cs.2018.5$y.2021 = 1
cs.2018.5$y.2022 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2018[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2018
CS.ng$y.2019[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2019
CS.ng$y.2020[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2020
CS.ng$y.2021[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2021
CS.ng$y.2022[match(cs.2018.5$CSREF,  CS.ng$CSREF)] <- cs.2018.5$y.2022
}



#> 4. Agreements starting in 2018 running for 10 years 
cs.2018.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2018-01-01") & STARTDATE.2 <= as.Date("2018-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2018.10)[1] != 0) {
#> Update relevant cols
cs.2018.10$y.2018 = 1
cs.2018.10$y.2019 = 1
cs.2018.10$y.2020 = 1
cs.2018.10$y.2021 = 1
cs.2018.10$y.2022 = 1
cs.2018.10$y.2023 = 1
cs.2018.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2018[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2018
CS.ng$y.2019[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2019
CS.ng$y.2020[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2020
CS.ng$y.2021[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2021
CS.ng$y.2022[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2022
CS.ng$y.2023[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2023
CS.ng$y.2024[match(cs.2018.10$CSREF,  CS.ng$CSREF)] <- cs.2018.10$y.2024
}



#> 5. Agreements starting in 2019 running for 5 years 
cs.2019.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2019-01-01") & STARTDATE.2 <= as.Date("2019-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2019.5)[1] != 0) {
#> Update relevant cols
cs.2019.5$y.2019 = 1
cs.2019.5$y.2020 = 1
cs.2019.5$y.2021 = 1
cs.2019.5$y.2022 = 1
cs.2019.5$y.2023 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2019[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2019
CS.ng$y.2020[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2020
CS.ng$y.2021[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2021
CS.ng$y.2022[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2022
CS.ng$y.2023[match(cs.2019.5$CSREF,  CS.ng$CSREF)] <- cs.2019.5$y.2023
}



#> 6. Agreements starting in 2019 running for 10 years 
cs.2019.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2019-01-01") & STARTDATE.2 <= as.Date("2019-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2019.10)[1] != 0) {
#> Update relevant cols
cs.2019.10$y.2019 = 1
cs.2019.10$y.2020 = 1
cs.2019.10$y.2021 = 1
cs.2019.10$y.2022 = 1
cs.2019.10$y.2023 = 1
cs.2019.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2019[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2019
CS.ng$y.2020[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2020
CS.ng$y.2021[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2021
CS.ng$y.2022[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2022
CS.ng$y.2023[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2023
CS.ng$y.2024[match(cs.2019.10$CSREF,  CS.ng$CSREF)] <- cs.2019.10$y.2024
}



#> 7. Agreements starting in 20 running for 2years 
cs.2020.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.2)[1] != 0) {
#> Update relevant cols
cs.2020.2$y.2020 = 1
cs.2020.2$y.2021 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.2$CSREF,  CS.ng$CSREF)] <- cs.2020.2$y.2020
CS.ng$y.2021[match(cs.2020.2$CSREF,  CS.ng$CSREF)] <- cs.2020.2$y.2021
}



#> 8. Agreements starting in 20 running for 5 years 
cs.2020.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.5)[1] != 0) {
#> Update relevant cols
cs.2020.5$y.2020 = 1
cs.2020.5$y.2021 = 1
cs.2020.5$y.2022 = 1
cs.2020.5$y.2023 = 1
cs.2020.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2020
CS.ng$y.2021[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2021
CS.ng$y.2022[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2022
CS.ng$y.2023[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2023
CS.ng$y.2024[match(cs.2020.5$CSREF,  CS.ng$CSREF)] <- cs.2020.5$y.2024
}



#> 9. Agreements starting in 20 running for 10 years 
cs.2020.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2020-01-01") & STARTDATE.2 <= as.Date("2020-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2020.10)[1] != 0) {
#> Update relevant cols
cs.2020.10$y.2020 = 1
cs.2020.10$y.2021 = 1
cs.2020.10$y.2022 = 1
cs.2020.10$y.2023 = 1
cs.2020.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2020[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2020
CS.ng$y.2021[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2021
CS.ng$y.2022[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2022
CS.ng$y.2023[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2023
CS.ng$y.2024[match(cs.2020.10$CSREF,  CS.ng$CSREF)] <- cs.2020.10$y.2024
}



#> 10. Agreements starting in 21 running for 2 years 
cs.2021.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.2)[1] != 0) {
#> Update relevant cols
cs.2021.2$y.2021 = 1
cs.2021.2$y.2022 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.2$CSREF,  CS.ng$CSREF)] <- cs.2021.2$y.2021
CS.ng$y.2022[match(cs.2021.2$CSREF,  CS.ng$CSREF)] <- cs.2021.2$y.2022
}


#> 11. Agreements starting in 21 running for 5 years 
cs.2021.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.5 )[1] != 0) {
#> Update relevant cols
cs.2021.5$y.2021 = 1
cs.2021.5$y.2022 = 1
cs.2021.5$y.2023 = 1
cs.2021.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2021
CS.ng$y.2022[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2022
CS.ng$y.2023[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2023
CS.ng$y.2024[match(cs.2021.5$CSREF,  CS.ng$CSREF)] <- cs.2021.5$y.2024
}



#> 11. Agreements starting in 21 running for 10 years 
cs.2021.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.10)[1] != 0) {
#> Update relevant cols
cs.2021.10$y.2021 = 1
cs.2021.10$y.2022 = 1
cs.2021.10$y.2023 = 1
cs.2021.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2021
CS.ng$y.2022[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2022
CS.ng$y.2023[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2023
CS.ng$y.2024[match(cs.2021.10$CSREF,  CS.ng$CSREF)] <- cs.2021.10$y.2024
}



#> 12. Agreements starting in 21 running for 20 years 
cs.2021.20 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2021-01-01") & STARTDATE.2 <= as.Date("2021-12-31") & Years == 20)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2021.20 )[1] != 0) {
#> Update relevant cols
cs.2021.20$y.2021 = 1
cs.2021.20$y.2022 = 1
cs.2021.20$y.2023 = 1
cs.2021.20$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2021[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2021
CS.ng$y.2022[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2022
CS.ng$y.2023[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2023
CS.ng$y.2024[match(cs.2021.20$CSREF,  CS.ng$CSREF)] <- cs.2021.20$y.2024
}



#> 13. Agreements starting in 22 running for 2 years 
cs.2022.2 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 2)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.2)[1] != 0) {
#> Update relevant cols
cs.2022.2$y.2022 = 1
cs.2022.2$y.2023 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.2$CSREF,  CS.ng$CSREF)] <- cs.2022.2$y.2022
CS.ng$y.2023[match(cs.2022.2$CSREF,  CS.ng$CSREF)] <- cs.2022.2$y.2023
}



#> 14. Agreements starting in 22 running for 5 years 
cs.2022.5 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 5)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.5)[1] != 0) {
#> Update relevant cols
cs.2022.5$y.2022 = 1
cs.2022.5$y.2023 = 1
cs.2022.5$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2022
CS.ng$y.2023[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2023
CS.ng$y.2024[match(cs.2022.5$CSREF,  CS.ng$CSREF)] <- cs.2022.5$y.2024
}



#> 15. Agreements starting in 22 running for 10 years 
cs.2022.10 <- CS.ng %>%
  select(CSREF, STARTDATE.2, Years, y.2016, y.2017, y.2018, y.2019, y.2020, y.2021, y.2022, y.2023, y.2024) %>% 
  filter(STARTDATE.2 >= as.Date("2022-01-01") & STARTDATE.2 <= as.Date("2022-12-31") & Years == 10)
#> If data frame has no observations skip the subsequent updating of main CS.ng table columns
if (dim(cs.2022.10)[1] != 0) {
#> Update relevant cols
cs.2022.10$y.2022 = 1
cs.2022.10$y.2023 = 1
cs.2022.10$y.2024 = 1
# Update main CS.ng table with data from agreement subset group
CS.ng$y.2022[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2022
CS.ng$y.2023[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2023
CS.ng$y.2024[match(cs.2022.10$CSREF,  CS.ng$CSREF)] <- cs.2022.10$y.2024
}



# Remove records where average annual cost id "unavailable"
CS.final <- CS.ng[!(CS.ng$AVGANNCOST=="*Unavailable"),]



# Export final raw data table
write_csv(CS.final, here("Out", "AES", "CS", "Nat_Park", "CSV", paste0("CS_", active_park, "_Raw_Processed.csv")))

# # Export final raw data table - TESTING
# write_csv(CS.final, here("Out", "AES", "CS", "County", "CSV", paste0("CS__Raw_Processed_Test.csv")))



# Change annual cost column to numerical
CS.final$AVGANNCOST <- as.numeric(CS.final$AVGANNCOST)
CS.final$AVGANNCOST <- round(CS.final$AVGANNCOST, 2)


#> Calculate annual costs
#> Total for 2018
CS.2018.Total <- CS.final %>% 
  filter(y.2018 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2018 <- sum(CS.2018.Total$AVGANNCOST)

#> Total for 2019
CS.2019.Total <- CS.final %>% 
  filter(y.2019 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2019 <- sum(CS.2019.Total$AVGANNCOST)


#> Total for 2020
CS.2020.Total <- CS.final %>% 
  filter(y.2020 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2020 <- sum(CS.2020.Total$AVGANNCOST)


#> Total for 2021
CS.2021.Total <- CS.final %>% 
  filter(y.2021 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2021 <- sum(CS.2021.Total$AVGANNCOST)

#> Total for 2022
CS.2022.Total <- CS.final %>% 
  filter(y.2022 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2022 <- sum(CS.2022.Total$AVGANNCOST)

#> Total for 2023
CS.2023.Total <- CS.final %>% 
  filter(y.2023 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2023 <- sum(CS.2023.Total$AVGANNCOST)

#> Total for 2024
CS.2024.Total <- CS.final %>% 
  filter(y.2024 == 1) %>% 
  select(AVGANNCOST)
#> Variable to hold total sumemd value for year
CS.sum.2024 <- sum(CS.2024.Total$AVGANNCOST)


#> Create data frame of annual costs
CS_Summary <- as.matrix(data.frame(Total_2018 = numeric(0), Total_2019 = numeric(0), Total_2020 = numeric(0), Total_2021 = numeric(0), Total_2022 = numeric(0), Total_2023 = numeric(0), Total_2024 = numeric(0)))
totals.cs <- c(CS.sum.2018, CS.sum.2019, CS.sum.2020, CS.sum.2021, CS.sum.2022, CS.sum.2023, CS.sum.2024)
#> Insert figure for 2018
df2 <- miscTools::insertRow(CS_Summary, 1, totals.cs)
df2 <- as.data.frame(df2)

#> Print table district.summary on web page
print(knitr::kable(df2, caption = paste0("Table 4.6.", i , " ", "CS payments - ",  active_park)))

#> Make long
df2.plot <- melt(df2)
#> Extract year from string
names(df2.plot)[1]<-paste("Year")
names(df2.plot)[2]<-paste("Value_cs")
df2.plot$Year <- sub("^.*([0-9]{4}).*", "\\1", df2.plot$Year)

# Simple bar chart
plot.cs <- ggplot(df2.plot, aes(x=Year, y=Value_cs)) + 
  geom_bar(stat = "identity", fill="#e38d0b", width = 0.6) +
  theme_bw() +
  labs(title = paste0("Countryside Stewardship payments: ", active_park), x = "Year", y = "Value (£)") +
  geom_text(aes(label= paste0(round(Value_cs / 1000000, digits = 2), " M")), vjust = 1.5, colour = "white", size = 3.5) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 16, b = 0, l = 0))) +
  theme(axis.title.x = element_text(margin = margin(t = 12, r = 0, b = 0, l = 0))) +
  theme(legend.position="none") +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))
plot.cs
print(plot.cs)


#> Export the plot
ggsave(here("Out", "AES", "CS", "Nat_Park", "Plot", paste0("CS_payments_", active_park, ".png")), width = 7, height = 5)

# #> Export the plot - TESTING
# ggsave(here("Out", "AES", "CS", "County", "Plot", paste0("CS_payments_County_Test.png")), width = 7, height = 5)


#Export table
write_csv(df2, here("Out", "AES", "CS", "Nat_Park", "CSV", paste0("CS_payments_", active_park, ".csv")))


# #Export table - TESTING
# write_csv(df2, here("Out", "AES", "CS", "County", "CSV", "CS_payments_County_Test.csv"))

}
Table 4.6.1 CS payments - Exmoor
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
201416.2 979555.7 1701834 2732174 2624228 2210807 1662666

Table 4.6.2 CS payments - Dartmoor
Total_2018 Total_2019 Total_2020 Total_2021 Total_2022 Total_2023 Total_2024
217712.2 530082.7 1180132 1988774 1852456 1484083 1171713

5. Agri land proportions calcs

Summary of Corine land cover types (all) by:

  • County level

  • District level

  • NPs

# #> Get corine data (pre-clipped to South West in QGIS
# corine.southWest <- st_read(here("In", "Shape", "Corine4_landStats.shp"))
# 
# #> Get land class description from lookup
# #> Get LUT
# corine.LUT <- read_csv(here("In", "Shape", "clc_legend.csv"))
# corine.southWest <- merge(corine.southWest